Accuracy of numerical wave model results: Application to the Atlantic coasts of Europe

. Numerical wave models are generally less accurate in the coastal ocean than offshore. It is generally suspected that a number of factors specific to coastal environments can be blamed for these larger model errors: complex shoreline and topography, relatively short fetches, combination of remote swells and local wind seas, less accurate wind fields, presence of strong currents, bottom friction, etc. These factors generally have strong local variations, making it all the more difficult to adapt a particular model setup from one area to another. Here we investigate a wide range of modelling choices including 5 forcing fields, spectral resolution and parameterizations of physical processes in a regional model that covers most of the Atlantic and North Sea coasts. The effects of these choices on the model results are analyzed with buoy spectral data and wave parameters’ time series. Additionally, satellite altimeter data is employed to provide a more complete performance assessment of the modelled wave heights as function of the distance to the coast and to identify areas where wave propagation is influenced by bottom friction. We show that the accurate propagation of waves from offshore is probably the most important factor on 10 exposed shorelines, while other specific effects can be important locally, including winds, currents and bottom friction.

exception of those placed less than 800 m from the coastlines, were included in the generation of the new mesh fixing their previous position. This was done to facilitate the use of the new results by users of the previous hindcast.
Finally, the resolution was increased in 14 zones of interest for marine energy users, (Fig. 1a). The generated mesh has a total of 328030 nodes (Fig. 1b), with a resolution (triangle side) ranging from ∼200 m at the coast and refined zones to, 60 approximately 15 km in deep offshore areas. An alternative to this careful editing of the mesh is the use of implicit schemes. However, using implicit schemes with CFL values much larger than 1 opens the door to both larger advection errors (stability does not imply accuracy) and larger splitting errors as the time steps for advection can be much smaller than the refraction and source term time step (Roland and Ardhuin, 2014). We have preferred to stick to the explicit N-scheme because numerical efficiency is not central in the study,

Bottom sediment map
The construction of a sediment grain size map was included to properly represent wave energy dissipation due to bottom 70 friction (see section 5.5 for results). In the model, the grain size is characterized by its median diameter D 50 , defined at each node of the mesh. The D 50 values where estimated from the EMODnet harmonized seabed substrate charts. The minimum grain size was set to 0.02 mm, while zones characterized as pebbles or larger elements (boulders) were represented with a D 50 =150 mm. By default, the minimum grain size was applied to all regions where no substrate was specified. Since most areas with no bottom characterization are in deep waters (e.g. > 400 m), this assumption does not have any relevant effect on 75 the wave fields evolution. The bottom sediment diameter map is presented in Fig. 2.

Source terms and numerical choices
In WW3 the wave action equation is solved using a splitting method to treat in different steps temporal depth changes, spatial propagation, intra spectral propagation and source terms (Yanenko, 1971;Tolman and Booij, 1998; The WAVEWATCH III ® Development Group, 2019). Spectral propagation, which includes refraction, is computed with an explicit third order scheme that combines the QUICKEST scheme with the ULTIMATE total variance diminishing limiter (Leonard, 1991), while spatial advection is done with the explicit Narrow stencil scheme (N-scheme) (Csík et al., 2002;Roland and Ardhuin, 2014). Non-linear evolution and wave to wave interactions are represented with the Discrete Interaction Approximation (DIA, Hasselmann et al., 1985). The utilized wind input and wave dissipation source terms are taken from the ST4 parameterizations described in Ardhuin et al. (2010) with adjustments described in Alday et al. (2021) consistent with the global model used for our boundary 85 conditions. A constant wave energy reflection of 5% is used at the coastlines, as parameterized by .
In the present study we only analyze the effects of changes in the ST4 parameterizations. A detailed list of the parameters used for the model implementation is given in Appendix A.

Boundary conditions and forcing fields improvements
The accuracy of modelled wave data directly depends on the quality of the forcing fields and the provided boundary conditions 90 (BC) for the case of nested models. This becomes particularly relevant in coastal areas, for accounting wave-current interactions in macro tidal areas, the assessment of energy resources, port design and operation conditions, or the study of extreme events.
Along with the high spatial resolution, an important aspect of the wave hindcast analyzed in this study, is the utilization of improved spectral BC from the wave data set described in Alday et al. (2021). This wave hincast was created using wind fields from the fifth generation ECMWF atmospheric reanalyses of the global climate, ERA5 (Hersbach et al., 2020), and surface 95 current fields taken from the CMEMS Global Ocean Multi Observation Products (MULTIOBS_GLO_PHY_REP_015_004).
The global grid from where boundary conditions are taken has a spatial resolution of 0.5 • , while the wave spectrum is discretized in 24 directions (15 • resolution) and 36 exponentially spaced frequencies from 0.034 to 0.95 Hz with a 1.1 increment factor from one frequency to the next. The proposed spectral discretization, wave growth and dissipation parameters, along with the use of upgraded forcing fields, showed clear improvements of sea state parameters (at global scale) when compared 100 to previous hindasts, like the widely used data set from Rascle and Ardhuin (2013).
The ( For the proposed regional model, three main forcing fields were included: wind, tidal levels and tidal currents. As for the 105 global model, ERA5 surface winds were used for wave generation. Similar to what was done in Boudière et al. (2013), tidal levels and currents time series were reconstructed in WW3 with harmonics taken, in this case, from two different sources. The first one, is the output from Ifremer's tidal atlas (Pineau-Guillou, 2013) created with MARS 2D (Lazure and Dumas, 2008), a hydrodynamic model based on the shallow water equations. A total of 5 embedded models with 3 levels of nesting and different spatial resolutions were selected (Fig. 3a). The second tidal data source was used to cover part of the Atlantic coast of Portugal 110 until the Gulf of Cadiz, which are not included in the tidal atlas. The complement data was taken from the native mesh of the FES2014 model (Carrere et al., 2015) and regridded to 0.004° (Fig. 3b).
In all simulations, the boundary conditions are updated every 3 hours, winds every 1 hour, tidal levels and velocities fields are updated each 30 minutes. The output frequency of the nested model is hourly.

115
The same extended frequency range used in the global grid was employed in the regional mesh to perform all simulations, matching the discretization at the boundary. The extension to higher frequencies is aimed to allow for a better representation of the variability of the wave spectrum for very low wind speeds or very short fetches. At the other end, the purpose of adding lower frequencies is to let the spectrum develop longer wave components for severe storm cases (e.g. Hanafin et al., 2012).
In terms of directional discretization, we used 36 mainly directions (10 • resolution), and tests with 24 and 48 directions were 120 employed to verify the effects of the directional resolution.
The source terms are integrated with an adaptative time step that is automatically adjusted in the range 5 to 180 s. We defined the maximum model advection time step to be 30 s, taking into account the minimum mesh triangle area and the presence of strong currents. The refraction time step was set to 15 s. Sensitivity tests with smaller values (not shown) had very limited impact on the model results.

Buoy data
We use 6 French buoys with spectral data provided by CEREMA, and 2 Belgian buoys from which spectra were not available, but besides the usual significant wave height, they provide a low frequency wave height H 10 (Fig. 4). The H 10 pararameter corresponds to a wave height computed for periods from 10 s and longer (≤ 0.1 Hz). These sites cover a wide range of depths, 130 current intensities, tidal amplitude levels and proximity to shore, which makes them an appropriate sampling group to evaluate the overall accuracy of the results (table 1). No assessment of potential instruments' replacements, maintenance periods nor deploy position changes have been taken into account for this study.
To match the frequencies discretization of the spectrum and output frequency (hourly) in WW3, spectral data from the in situ measurements have been first interpolated into the same discrete frequencies used in the model, and then averaged in time 135 to provide hourly output.

Buoy
where X obs are the observed quantities from in situ or satellite measurements, X mod are the modelled quantities (spectral values or integrated wave parameters), and N the total amount of analyzed data.
We use the term normalized mean differences (NMD) when using eq. 5 between different model configurations.
5 Sensitivity analyses results and discussion 155

Influence of spatial resolution
We first consider significant wave heights (hereinafter H s or "wave height"). A comparison between February 2011 mean H s fields from the global model described in section 2.4 and our implemented regional model is presented in Fig. 5. To evaluate the differences between models, the output from the 0.5 • global grid was linearly interpolated onto the regional mesh nodes before computing the mean wave height and the mean difference for the selected time window.

160
The most important differences are found on the shelf where complex coastline geometry and bathymetry requires higher detail to better represent land shadows and wave refraction (NMD in Fig. 5). The largest positive differences (> 20%) are commonly found in the regions sheltered from North Atlantic swells. In the global model, islands and headlands smaller than the grid size are parameterized as obstructions of the wave energy flux (Chawla and Tolman, 2008). Another direct effect of using increased spatial resolution can be seen between the Orkney and the Shetland islands. The regional model shows 165 averaged wave height values of almost 5 m in this area for the analyzed month. On the other hand, the combined effects of the sub grid obstruction approach and coarse resolution of the global grid, leads to high underestimation of about -20% with respect to our mesh results.

Regional model mean Hs
Hs NMD (Global -Regional) Figure 5. Mean wave height fields from global and regional model, and wave height normalized mean differences (NMD: Global -Regional).
Dashed black lines represent the 400 m depth contours. Areas where no wave data are available from the global grid are highlighted with a gray background in left and right panels. Results for February 2011.

170
(2010) and Leckler et al. (2013), these adjustments were designed to better represent the wave heights measured by altimeters at global scales. Here we further consider the impact of these modifications on waves in our coastal domain, using five different simulations with parameter changes listed in table 2. These changes include an empirical enhancement of the wind speeds above a threshold U c by the amount x c (U 10 − U c ), and a modification of the swell dissipation with a change in the threshold Reynolds number Re c that defines the transition from the weak (laminar) to strong (turbulent) swell dissipation and the swell 175 dissipation coefficient s 7 .
We analyzed model results for two months when extreme sea states have been recorded, February 2011 and January 2014. In February 2011, the extra tropical storm Quirin generated extreme sea states with peak periods exceeding 20 s over the western coasts of Europe. In January 2014, storm Hercules was one of the many storms from a particularly severe winter. This event caused vast coastal damage in the United Kingdom (Masselink et al., 2016), and from the Western coast of France to Portugal 180 (Masselink et al., 2015). Wave height values exceeded 10 m and peak periods exceeded 20 s (Ponce De León and Soares, 2015; Castelle et al., 2015). Given the characteristics of the selected cases, it is considered that they are suitable to study wave energy Despite the similarities between time series of the wave parameters such as H s and T m02 from one test to another, they 185 noticeably differ for extreme values. Yet, the model runs have systematic differences as a function of wave heights or wave periods, with 5 to 10% deviations for the larger periods and heights that correspond to the most severe storms and associated swells ( Fig. 6). In these events, and consistent with the global scale results, the wind enhancement is most effective at correcting the low bias in extreme wave heights and mean periods that is typical of the previous hindcasts. Adjustments to the swell dissipation have a negligible impact when acting only 1000 km or less of propagation within our coastal domain. As shown in 190 Fig. 7, the wind enhancement allows the generation of lower frequency waves. This improves the model accuracy at exposed buoys 62066, 62074 and 62069, and produces realistic energy levels for frequencies below 0.05 Hz during the extreme events of January 2014. Unfortunately, the correction also produces too much low frequency energy at the shallower buoy 62078. We suspect that dissipative processes in shallow water may be underestimated for these very large periods (Fig. 7e,f).   195 At global scale, the use of ocean surface currents can improve the accuracy of the simulated sea states (Echevarria et al., 2021;Alday et al., 2021), although a full effect generally requires relatively high spatial resolution that is generally not achievable by observations and thus models are usually not constrained at the necessary scale (Marechal and Ardhuin, 2020). This is the main reason why geostrophic currents were not considered in the high resolution regional model.

Wave-Current Interactions
Adding surface currents in the simulations can have effects on wave generation due to changes on the relative wind, it can 200 modify the advection of waves or induce refraction in regions with large current gradients. Given the diverse tidal amplitudes within the modeled domain, it is expected to have different effect levels over the sea states in different areas. We thus attempt to characterize the changes of the wave field when tidal currents are taken into account in the simulations. To do so, we look at differences on a set of wave parameters, namely H s , directional spreading SPR, the peak direction D p and peak period T p .
We first checked global scale current effects via the boundary conditions, and then focused on tidal current effects within our 205 coastal domain.
To evaluate the effects of global currents on the boundary condition, we analyzed a specific output time with a large Atlantic swell, and differences between 1 month simulations. The most noticeable changes caused by global currents are obtained for H s , D p and directional spreading ( Fig. 8 middle panel), with typical differences of the order of 5 %. These differences vanish when averaged over one month (Fig. 8 right panel).  The effects of tidal currents within the model domain are generally more important, with some strong local effects caused by the high spatial currents' variability. In contrast to the influence of global currents in the BC, there is a clear increase of the wave fields' differences at each temporal output, that can be larger than +/-10 %. Feature mainly seen along the English Channel and the Irish Sea (Fig. 9, left panel). Over the entire month, tidal currents induce mean wave height differences of the order of 5 % (Fig. 9, right panel). 215 The use of tidal currents also proved to have large impact over the peak period (T p ), up to 15% differences in Normandy and Liverpool bay, for example, and 8 % mean differences over one month (not shown).
There is a noticeable feature of the wave field along the shelf break, starting at the Bay of Biscay and extending northwards up to 49°N, which can be seen more clearly through the D p and H s fields from Fig. 8b,c (left panel), and particularly by analyzing the effects of tidal currents over the wave heights in Fig. 9c (left panel). The intensities of current used in our model 220 present maximum values of about 0.5 m s −1 along the aforementioned area, which is consistent with previously recorded in situ measurements and the expected sharp variation of currents across the shelf break (i.e. Le Cann, 1990). It is thought that the distinct gradients visible in some of the wave parameters are function of the tides' phase and the mean wave direction.
Attempts to identify the presence of this signature with altimeter data is an ongoing subject of study.
Results were further compared against in situ data from January 2014 at buoy 62059 (Fig. 10). Including tidal currents helps 225 to reduce the high energy bias at low frequencies, probably due to an overall reduction of the effective wind input for locally generated waves during the tidal cycle (Fig. 10a). In Fig. 10c is possible to observe the modulation of H s and T m01 caused by the changes in currents intensities and direction (blue line in figure), which in the end helps to reduce the bias of these quantities compared to the measurements (Fig. 10b). Notice that there is a constant shift in the occurrence of peaks and troughs of H s and T m01 in Fig. 10c. This is thought to be mostly attributed to a slight phase shift in the tidal forcing field, which introduces

Effects of spectral directional resolution
The selection of the spectral discretization plays an important role in the characteristics of the modeled sea states (Tolman, 1995b;Roland and Ardhuin, 2014). Normally, in coastal applications like assessments of wave energy or simulation of storm surges, higher time and spatial variability details are desired, and hence, higher spatial and spectral resolution is required (e.g. 235 Bertin et al., 2015;Accensi et al., 2021;Wu et al., 2020). Nevertheless, the quality of the results may be affected by the characteristics of the used BC.
We analyzed the changes in the energy distribution of the directional spectrum and the wave field evolution due to different directional resolution values in our mesh and in the BC. The different BC tests are aimed to identify potential effects when coarser resolution is used at global scale, and then interpolation is applied to match the resolution of the nested mesh ( Table 3. Tests for spectral directional resolution effects. All parameters not specified here correspond to test T475. When directional resolution of the boundary conditions (BC) is lower than in the mesh, interpolation is applied at the boundary to match the resolution of the nested model.
Variations in the energy distribution due to lower resolution in the BC are presented in Fig. 11, comparing BC with 24 spectral directions with respect to 36. A set of 4 locations were selected: At the boundary (named node W12N56), and along 250 the French coast nodes 62074 (Belle Ile), 62069 (Pierres Noires) and 62059 (Cherbourg). Bathymetry details of these locations were presented in Fig. 4. Here we present results for January 2014, but the analysis was also conducted with simulations for February 2011.
At the boundary, most of the difference in energy traveling outside the domain are related to very low levels of spectral energy (angles > 270°and < 360°, Fig. 11a right panels). This has negligible effects over of the already analyzed wave parameters 255 (e.g. H s , D m , SPR). For waves traveling into the domain , only large differences (N M D > +/ − 10 %) are observed at lower frequencies (< 0.1 Hz) between directions 20°and 150° (Fig. 11a right panel), which corresponds to the area with higher mean energy at this location for January 2014 (defined by the contours in Mean Energy panel of Fig. 11a). We found that this effect is still present in nearshore areas exposed to the incoming swells from the North Atlantic (nodes 62074 and 62069), although with an overall narrower directional range attributed mostly to the bathymetry induced refraction that tends to "align" the arriving 260 waves (Fig. 11b,c).
No significant changes in energy distribution were found at node 62059, for each output time and for the full simulation ( Fig. 11d). This is expected since at Cherbourg the sea state characteristics are mostly driven by the local winds.
To further assess potential changes introduced in wave parameters, we analysed the differences in fields of H s , T p , SPR, D p , and the mean direction D m (not shown; Fig. 12). Using coarser directional resolution in the BC has minor effects over wave 265 parameters integrated along the complete frequency range (e.g. D m or H s ; Fig. 12b, top panel). Differences in the results are exacerbated when BC with 24 directions are interpolated into 48 (right panels in Fig. 12a,b) but in general mean and random differences between tests remained below +/−2.5 %, with the exception of T p that presented the largest NRMSD.
Even though the magnitude of these quantities remain fairly consistent, interpolating BC with coarser directional resolution affects the characteristics of the wave fields propagating into the domain. This is attributed to slight changes in the wave 270 celerity (C=gT /2π in deep waters) due to frequency shifts in the neighborhood of the energy peak ( Fig. 11a,b,   The analysis of directional resolution of the mesh is mainly focused on the Garden Sprinkler Effect (GSE) on wave propagation. This phenomenon is observed as a separation or disintegration of continuous swell fields propagated with a discretized spectral wave model (Booij and Holthuijsen, 1987;Tolman, 2002). The generation of the GSE is namely linked to the spectral 275 resolution and the advection scheme. Currently there are no GSE alleviation methods available for unstructured grids in WW3.
A good example was found during February 1st 2011, where a strong swell from the North Atlantic arrived to the northern coast of Scotland. In Fig. 13a we present an instant (13:00:00 UTC) of the event using 3 different discretizations from tests 24D24BC, 36D36BC, 48D48BC (table 3). The GSE can be observed to the East of the Orkney and Shetland Islands towards the Norwegian Sea (between latitudes 58°and 61°) when 24 directions are employed (Fig. 13a, left panel).

280
The impact of the GSE was assessed by comparing the results against the output from a model with higher directional resolution. Via a straight forward difference between tests, is possible to visualize changes of the H s field caused by the spurious wave propagation pattern (Fig. 13b). Comparing tests 24D24BC with 36D36BC, and for this particular scenario, differences in wave height can reach +/-0.2 m (roughly +/-5%) as the swell approaches Norway, between longitudes 2°to 4°( The mild repercussion of the GSE over the wave height field in the present results shouldn't be generalized, since this phenomenon could be intensified depending on the incoming swell conditions. Our findings suggest that using a directional resolution of 10 • is enough to mitigate the effects of the GSE. It is relevant to point out that, for example, the required computation time in 36D36BC is 40% higher than in 24D24BC, a considerably elevated cost for potential operational (forecasting)

Bottom friction effects
Over the continental shelf, in intermediate to shallow waters, the evolution of the wave fields becomes influenced by the bottom characteristics. In the absence of strong wind seas and outside the surf zone, dissipation of energy is mainly induced by bottom roughness effects. We thus try to quantify the effects of including the bottom friction sink term in the wave action equation.

295
To provide a general view, we compared model output from 1 year tests with the 1 Hz altimeter data from the ESA Sea State CCI V2 datset. For this particular analysis 1 year simulations were required in order to have at least a minimum of 5 satellite measurements to compare with the regridded WW3 wave height fields at 0.1°. Only altimeter measurements at least 10 km away from the coastline were considered to avoid potential data with high noise to signal ratio.
Bottom friction effects were included through the SHOWEX parameterization proposed in Ardhuin et al. (2003). This 300 expression was initially developed for sandy bottoms based on the eddy viscosity model of Grant and Madsen (1979) and includes a decomposed roughness parameterization for ripple formation and sheet flow. In WW3 it has been implemented including a sub grid parameterization for water depth variability following Tolman (1995a). The bottom friction source term can be written as follows: and When the Shields number ψ is >= A 3 ψ c , the Nikuradse roughness k N is taken as the sum of the ripple roughness k r and a 310 sheet flow roughness k s : where, with, In table 4 we present a set of empirical parameters originally taken from Ardhuin et al. (2003) where we have particularly modified A 5 to 0.04. D 50 is the median sediment size in meters defined at each node of the unstructured grid (see Fig. 2). A full description of the terms in eq. 6 to 13 can be found in Ardhuin et al. (2003)  To assess the effects of the bottom friction parameterization, we first compared 1 year simulations with and without dissipa-325 tion to verify changes in the wave field. In Fig. 14a we present the wave height mean bias obtained by comparing with Saral (year 2014) for the full domain. A clear reduction of the wave heights bias is detected in the south of the North Sea. In this area, we found that wave height mean differences between results with and without bottom friction can be of 0.3 m and higher.
Analysis with other altimeters (e.g. Jason-2 and Envisat) for year 2011 show consistent results.
In general, with altimeter data most relevant changes in wave heights, when bottom friction is included, are detected for 330 depths smaller than 50 m. We found a couple of Envisat tracks passing almost parallel off the coast of La Rochelle and close to Ile de Yeu (Fig. 14b). In both locations the use of the bottom friction parameterization, with the defined D 50 , helps to reduce the mean bias on wave heights. These results are consistent with the findings of Roland and Ardhuin (2014) for this area based on buoy data.
We picked 3 locations to compare our results with in situ measurements, buoy 62078 on the Altantic French coast, and buoys 335 Westhinder and Scheur Weilingen deployed in shallower depths along the coast of Belgium (Fig. 4).
For buoy 62068 we first compared the full time series of in situ wave height against simulations with and without bottom friction effects. Reductions in the wave height bias and NRMSE of respectively 4.5 % and 5.0% are obtained when bottom friction and the sediment size map are included (Fig. 15a,c). Nevertheless, most significant changes in the modeled wave height appear at wave heights roughly larger than 3 m. We then selected an ad hoc wave height threshold of 3.5 m to define 340 "extreme" sea states and analyze the effects of the parameterization over the events on which dissipation due to wave-bottom interactions is dominant. For these events, a wave height bias and RMSE reduction of about 0.3 m, with a decrease of about 8% and 5.3% in the bias and scatter is obtained when the SHOWEX dissipation term is used (Fig. 15b). Moreover, we found good agreement between the occurrences of the Shields number ψ exceeding its critical value ψ c (Fig. 15d) and the occurrences of extreme sea states with H s > 3.5 m (Fig. 15c) especially between January and March 2014. In the model, the evolution time 345 scale due to bottom friction is inversely proportional to f e u b,rms , which gives a measure of the strength of bottom friction, sharply increasing every time the critical Shields number is exceeded. In this case, the definition of extreme events helps to identify when the effects of bottom friction becomes relevant, since larger H s are normally related to longer wave lengths, thus wave-bottom interactions start at deeper depths.
At Westhinder and Scheur Weilingen we analyzed the dissipation effects over components of the spectrum with periods 350 longer than 10 s comparing H 10 values. For these locations we also compare with simulations using the JONSWAP bottom friction parameterization (Hasselmann et al., 1973;Tolman, 1991) with its default values (Fig.16). Wave energy for components longer than 10 s is clearly over estimated when no bottom friction is taken into account. The effect is visible at both analyzed depths. At Westhinder both parameterizations have similar effects, but at the shallower buoy location (Scheur Weilingen) the use of SHOWEX and the selected D 50 introduce a negative bias of H 10 > 0.5 m which could be related to an overestimation 355 of the sediment mean size in this area. WW3's frequency spectrum following eq. 6 to 13. D50 taken from bottom sediment map (Fig. 2)  Satellite altimetry provides a unique resource of worldwide wave heights' measurements. The integration and inter calibration of past and ongoing missions have allowed to continuously extend the coverage of measured data in space and time (Ribal and Young, 2019;Dodet et al., 2020). These datasets have been commonly used in open ocean applications to improve our 360 understanding of the sea states globally. On the other hand, their application in coastal (especially nearshore) areas has been very limited due to increased noise levels in the return signal. What is considered as noise is actually the detection of the non Gaussian land surface, which makes it difficult to retrieve the waves geophysical signal in the radar footprint.
The Sea State CCI V2 dataset employs the WHALES partial waveform retracking algorithm, more effective for reducing the intrinsic noise of the return signal, and suitable for coastal applications (Schlembach et al., 2020;Passaro et al., 2021). The 365 vast amount of measurements available at distances from the coast lines down to 5 km and less implies also a large coverage of measured wave heights in shallower depth areas, providing a broader description than local in situ records. Making use of the coverage and improvements in this altimeter product, we analyzed the performance of our mesh over part of the wave hindcast described in Accensi et al. (2021) which was created using the same mesh employed in the present study.
We analyzed 3 zones of the modeled area: Bay of Biscay, North Sea and English Channel. The purpose of the defined 370 zones is to assess the performance of the model in different wave generation and propagation conditions. The Bay of Biscay is constantly exposed to swells radiated from the North Atlantic. At the North Sea, wave conditions are dominated by the local winds blowing over a well defined fetch and partially influenced by the swells from the Norwegian Sea. Finally, at the English Channel, most of swells' energy arriving from the North Atlantic is blocked, refracted and dissipated on its western end, local waves are generated over a very short fetch, and it is highly influenced by its tidal regime. From distances to the coast of 15 km and more we noticed a constant positive bias ranging from 2 to 6% in the Bay of Biscay, and in some cases going up to ∼8% in the English Channel. At the North Sea bias changes are more constrained between +/-2% (Fig. 17a). The positive bias in the Bay of Biscay is thought to be related to the BC obtained from the global hindcast using T475, which was calibrated with the Jason-2 data from CCI V1. This data was indeed found to overestimate wave heights 385 recorded by offshore buoy measurements (Dodet et al., 2020), which has been corrected in V2. The English Channel stands out as high bias and scatter area which is thought to be caused by the reduced amount of valid altimeter measurements in this area and inaccuracies of the forcing fields. Finally, less influenced by the BC and with an extended fetch for wave growth, the North Sea presents the lowest bias values, which along with the lower scatter (Fig. 17b) shows the good performance of the proposed parameterization and model setup in this area.

390
An overall bias decrease is observed for distances to the coast smaller than 15 km, which implies that in general the wave heights from altimeters are higher than those from the modeled. Differences that are particularly more accentuated at the Bay of Biscay at distances from the coast shorter than 10 km. Even with the higher uncertainty of modeled/measured wave heights closer to the shore, the available altimetry data down to ∼6 km offshore still provides unprecedented access to coastal information that, even at this early stage, allows to evaluate the model performance. In the present study we investigated the drivers of model errors in coastal areas, and how choices of parameterization, forcing, spectral and spatial resolution, and boundary conditions affect the characteristics of the simulated sea states. Extensive sensitivity analyses were carried out with a high resolution regional wave model for European coastal waters using the WAVEWATCH III framework. The performed tests and analyses were aimed to assess when and where the choices in the model setup have 400 a significant effect in regions where wave interactions with complex bathymetry, tidal currents and bottom roughness become important in wave propagation.
Overall, spatial resolution is one of the most important elements in shallow depth areas. We found that higher spatial resolution adequate to solve bathymetry features and explicitly solving coastlines can introduce changes in modelled wave height of about 20% when compared to lower resolution global models. Differences become more significant below 400 m depth, in 405 areas where refraction and diffraction are dominant, or in regions sheltered from the most frequent swell conditions.
Changes in the energy distribution of the spectrum were analyzed mainly from two points of view, introduced by modifications in the parameterization, and due to changes in directional resolution. Modification of the swell dissipation terms did not impact significantly the wave energy distribution in the regional domain, although their effect become important at global scales . In general, the applied enhancement to intensities higher than 21 m s −1 in the ERA5 wind fields 410 improves the model accuracy at swell-exposed locations, helping to reproduce realistic energy levels for frequencies lower than 0.05 Hz, partially solving their otherwise high under estimation (more than -50% in some cases). These findings suggest that the considerations taken to generate the boundary conditions at global scale, are one of the most important factors on shorelines exposed to waves from the North Atlantic.
Spectral energy differences due to directional resolution choices are larger than 10% at frequencies lower than 0.1 Hz. Effect 415 that is visible from the boundary to the nearshore in zones influenced by the BC. Differences in wave parameters (SPR, T p , D p ) observed between model tests suggest that the proper selection of directions to define the BC and within the nested model will help to reduce random errors. It was also found that with 10°resolution, the GSE is successfully alleviated in the mesh.
Within areas with large tidal amplitudes, including tidal forcing (currents and levels) typically change wave parameters by about 10% at each output time, and locally much more (e.g. . These differences are reduced for H s and 420 D p for a monthly average, but can still be larger than 5% for the SPR and T p . These findings imply that even if the average wave heights might be well estimated without tidal forcing, the propagation and evolution of the wave fields will be different. This can be observed in the H s and T m01 time series at buoy 62059 (Fig. 10).
Comparing with wave heights retrieved from altimeter data with 1 year simulations, we identified areas influenced by bottom friction dissipation by looking at changes in H s . We found that these changes can be observed at depths smaller than 50 m. In 425 shallower areas of the North Sea and some sections of the Atlantic coast of France, including the SHOWEX bottom friction parameterization helps to reduce the H s bias. Comparisons between model and in situ measurements of H 1 0 revealed an underestimation of the wave energy in the low frequency bands in very shallow areas. This effect could be related to a higher sensitivity of the SHOWEX parameterization in very shallow depths, thus, dissipation induced in longer wave components is over estimated with our current model setup.

430
Using 5 available missions from the Sea State CCI V2 dataset we performed a validation of the modelled wave height as function of the distance to the coast, between years 2002 to 2018. We observed an overall increase of wave height differences with our model for distances to the coast smaller than 10 km that can reach -8% (in average) at 5 km from the coast. These differences are likely due to increased uncertainties in altimeter measurements within the last 10km from the coast, where coastal features are known to strongly impact radar waveforms (Vignudelli et al., 2019). 435 We found that in many cases time averaged differences between model setups or with respect to in situ data are small, but these differences can be significant at each output time, implying that the time evolution of the sea states is in fact different.
This could partially explain cases with low bias and still larger random errors (e.g. SI) in some locations, when modelled wave parameters are compared with measurements.
Due to the different characteristics of the modelled domain (e.g. bathymetry features, bottom sediment type, fetch and tidal 440 amplitudes) the factors driving the accuracy of the model cannot be completely generalized. Instead, through the proposed analyses we have identified where changes in the wave field characteristics are more significant with different choices in forcing, resolution and parameterizations. Yet, it is not straightforward to assess how the combination of these choices can potentially compensate errors in the simulations. We find that boundary condition effects are most easily evaluated at deep water or partially sheltered locations (see also Crosby et al., 2017), while separating bottom friction from other effects will 445 require a further analysis of specific swell events.
Data availability. The used coast line polygons, bathymetric data, bottom sediment type maps and buoy data have been take from the following web portals: -  Tolman, H. L. and Booij, N.: Modeling wind waves using wavenumber-direction spectra and a variable wavenumber grid, Global Atmos.