Current interaction in large-scale wave models with an application to Ireland

A nested large-scale wave model that covers most of the Atlantic Ocean and focuses on generating accurate swell conditions for Ireland is presented and validated over the two years 2016 and 2017. The impact of currents is studied using the surface currents from the GLORYS product. Currents slightly reduce by 1\% the error of the model compared to altimetry data, but they explain most of the wave energy at scales less than 50 meters. The variability induced by currents is found to be more noticeable on the instantaneous fields, up to 50 cm. Track following observations indicate that wave refraction induced by mesoscale eddies is correctly captured. However, this modulation is of the same order as a shorter-scale variability appearing in the altimetry data, which makes it difficult to objectively assess the impact of currents.


Introduction
Traditionally, large-scale simulations or hindcasts have been focusing on the air-sea interaction, trying to balance correctly wind input and swell dissipation. More recently, there has been some interest in the interaction with currents, which for the most part highlights the refraction induced by current gradients. The impact of large-scale currents like the Agulhas current or the Gulf Stream is well documented through numerical simulations cross-checked with satellite observations -see for example Marechal and Ardhuin, 2020;Barnes and Rautenbach, 2020 for recent studies on the Agulhas current. The impact of small scale currents has drawn less attention but in Ardhuin et al., 2017 it is shown that even submesoscale structures with shorter length scales (ranging from 10 km to 100 km) can explain the majority of the significant wave height spatial variability.
In wave spectra models, the interaction of currents on wave propagation appears in the wave action balance equation as a correction of the wave group velocity by a current velocity term. The framework used to derive this equation is usually a depth-uniform current, see for instance Peregrine, 1976. However, in the real world, the current is rarely uniform with depth. Therefore approximations are made, which differ depending on the application. A few are listed below.
Global wave models are computationally expensive. For this reason the current interaction is often ignored as the wind input and the bathymetry are the most important factors -see for example Stopa, 2018. In that paper, the wind is calibrated to the wave growth parameter within the spectral wave model WAVEWATCH III ® for 10 reanalysis datasets and 2 datasets composed of merged satellite observations. It is demonstrated that the space-time distributions of extreme waves are very different even after calibration. But wave-current interaction is ignored. The European Centre for Medium-Range Weather Forecasts (ECMWF) and the Copernicus Marine Environment Monitoring Service (CMEMS) each run a different forecast wave model based on the Wave Model numerical code (WAM). None takes currents into account. However both the ECMWF and CMEMS provide a wave reanalysis (ERA5 and WAVERYS respectively) coupled with an ocean model that takes into account the effects of currents on wave propagation. Large-scale wave models can have a more refined grid where the impact of current is more significant. In Marechal and Ardhuin, 2020, Ardhuin et al., 2017, and Barnes and Rautenbach, 2020 an offline coupling is done using an external current data-set. All these authors are directly interested in the wave height spatial variability caused by ocean currents. Different products are used for the currents but in each case they feature the surface current velocity. It is a cost efficient choice and still a good approximation in deep water (Kirby and Chen, 1989;Banihashemi et al., 2017). In a similar fashion the CMEMS runs a forecast model for the European North-West Atlantic shelf using WAVEWATCH III ® forced among others by surface ocean currents from one of their forecast ocean models.
The work presented in this paper consists of a large-scale application making use of an offline coupling to compute the current induced interaction terms. A standalone WAVEWATCH III ® model is set-up for the Atlantic and the resolution is refined around Ireland with two added nested grids. A sensitivity analysis on the wind input parameters is carried out first. It is shown that there is no single way to optimise the wind input growth term as the distribution does not react homogeneously to the growth parameter. The current interaction term in the wave action balance equation is computed. The current field corresponds to the surface currents from the 3D global ocean reanalysis GLORYS12 provided by the CMEMS. The impact of currents is evaluated by comparing the output of the model against satellite altimetry and wave buoy data, focusing on the North-East Atlantic region around Ireland. It is shown that currents can locally have a noticeable effect on wave propagation. However, they marginally impact the accuracy of the model on average. In Section 2, the wave model is described. Section 3 explores the sensitivity to wind input. The impact of the current field is investigated in Section 4. Conclusions are provided in Section 5.

Description of the wave model
A large-scale standalone WAVEWATCH III ® model (v6.07) is set-up and the different choices for the grid and model parameters are listed below. The wave action balance that is usually implemented in the most recent third generation wave models, like WAVEWATCH III ® , assumes a vertically uniform current that can slowly vary in the horizontal directions (the short-wave or large-scale current approximation). A complete derivation is given in Peregrine, 1976 and the same type of derivation is found in the more recent literature like Holthuijsen, 2010. The wave variables are denoted without a prime while the relative wave variables are denoted with a prime, except for the frequency where σ is the intrinsic frequency and ω the absolute frequency, related through the following dispersion relation highlighting a Doppler shift caused by the currents: (2.1) The following definitions are used for the absolute wave phase speed c and wave group velocity c g = (c g,x , c g,y ). Their intrinsic counterparts are derived by using the dispersion relation, with U = (U x , U y ) the current field passed down to the wave model and k = (k x , k y ) the wave vector with k its modulus: (2.2) The wave action balance is given by Eq. (2.4) below, with A(x, y, t, k, θ) the wave action, (x, y) the eastward and northward horizontal coordinates, t the time, σ the intrinsic wave frequency, θ the spectral direction, c θ and c k the rates of change of the direction θ and wavenumber k respectively, h the mean water depth, and S a global source term: (2.4) WAVEWATCH III ® (The WAVEWATCH III R ® Development Group, 2019) is a widely used third-generation wave model accounting for a wide range of processes. The source terms used in this paper are gathered in S (see Eq. (2.4)) and discussed briefly in the next paragraph, but a more exhaustive list is available in the WAVEWATCH III ® manual. It was shown to give excellent results for global simulations and is used in many operational models or reanalysis products (Ardhuin et al., 2011;Rascle and Ardhuin, 2013;Stopa, 2018). Different parameterizations are offered either for the propagation terms (left-hand side of Eq. (2.4)) or for the source terms (right-hand side of Eq. (2.4)). For the present application the source terms can be written as in Eq. (2.8) and some impactful choices are explicitly detailed in the next section: where the terms denote, respectively, linear growth, nonlinear transfer of wave energy through three-wave and four-wave interactions, wave growth by the wind, wave decay due to whitecapping, bottom friction and depth-induced wave breaking.

Current effects
The equations introduced above can only be derived for a vertically uniform current. Assuming a more realistic sheared current, these equations can still be used but the current field U becomes an effective current field, which in reality is a higher-order correction for the wave group velocity as explained in Chen, 1989 andBanihashemi et al., 2017. For this paper, the surface currents are used, which is common practice for large-scale wave models in deep water.
As shown in Eqs (2.4)-(2.1), currents are impacting wave propagation at different levels, and they can also impact the source terms in S. In Ardhuin et al., 2012, the main current impacts are shown to be the current-induced refraction through c θ , and the relative wind effect where the wind speed is taken relative to the currents in the wind input source term. It is worthwhile to note that similarly to the effect of current-induced refraction through c θ , the other current effects (current advection and Doppler shift) also contribute to increase the wave energy, and therefore wave heights, in areas of stronger opposing currents.

Grid set-up and parametrization
Three regular grids are used. They become more and more refined as one gets closer to Ireland, which is the area of interest for this study. The goal is to capture all the storms possibly impacting the wave conditions in Ireland, so the model does not require wave boundary conditions but only wind input. The grid specifications are summarized in table 1 with a graphical representation given in figure 1. The same common values are used for the three associated spectral grids and the spatial ratio between a child and a parent grid does not exceed 5, which is a commonly accepted value. The grid configuration is relatively similar to a previous wave hindcast focusing on the North-East Atlantic Ocean (Pilar et al., 2008).  Table 1: Grid specifications used in the WAVEWATCH III ® model. The goal is to obtain the best resolution possible for the North-East Atlantic grid, while conserving a similar number of grid points between all grids to balance the computational load.
The source terms that are included in this paper are highlighted in Eq. (2.8). For the most part they consist of standard choices suited for a large-scale application. The linear growth term S ln is parameterized as in Cavaleri and Rizzoli, 1981. It is useful to make the wave field grow from calm conditions. The nonlinear wave-wave interactions S nl are modeled with the Discrete Interaction Approximation from Hasselmann and Hasselmann, 1985. The source package ST4 (Ardhuin et al., 2010) is used. It encompasses the wind input S in and wave dissipation S ds . The impact of bed roughness appears through the bottom friction term S bot with the parametrisation of Grant and Madsen, 1979, and the depth-induced breaking S db is parameterized using Battjes and Janssen, 1978.

Wind input parametrization
A first sensitivity analysis is performed to adjust the wind input parametrization. Default values are provided in Ardhuin et al., 2010, but it is always recommended to adjust the most impactful parameters for a specific application (Ardhuin et al., 2011;Stopa, 2018). The choice of parameters for the wind input notably depends on the wind input source used and on the region for which the model is validated. The error on the output is not spatially homogeneous and different biases can be observed for different parts of the ocean (Ardhuin et al., 2011). These biases can be locally corrected by adjusting the model parameters. The wind input source term is part of the source package ST4 (Ardhuin et al., 2010) that also encompasses dissipation by whitecapping, which is an energy flux from the wave field to the ocean currents due to wave breaking not induced by the depth, and swell dissipation, which is a flux of energy from the wave field to the wind. Going quickly through the theory of wind-wave growth and referring here to the textbook Ardhuin, 2020, we emphasize that if the initial growth of waves is well explained by the turbulence theory developed in Phillips, 1957, the further growth of waves comes from the feedback of the waves on the atmosphere. The wave-induced pressure is roughly in quadrature with the surface elevation and that phase shift enables the transfer of energy from the air to the wave field at the surface. Using the linear equations from Airy's wave theory and adding an atmospheric term in quadrature with the surface elevation gives the simplified evolution of the wave amplitude a, shown below in Eq. (2.9), where β is a growth parameter. The right-hand side is then generalized for a wave system giving the spectral wind input source in Eq. (2.10), with N = σA the wave energy spectrum: (2.9) (2.10) The β term in Eq. (2.10) can be modeled by different formulations. Several options are available in WAVEWATCH III ® through the choice of the source package, and those parametrizations keep being refined as more and more field observations, experiments and high-resolution numerical simulations are made. The source package ST4 is taking into account several new features. For instance it is observed that the interaction between the air flow and the waves is stronger for younger waves, and at the same time the detachment of the air flow occurring for high winds is decreasing the wave growth, which can explain why the drag coefficient is observed to be reduced during hurricanes. A sheltering effect from the long waves on the short waves is also expected to reduce the growth of short waves. The wind input is computed as shown in Eq. (2.11) with ρ a and ρ w the air and water densities, κ the von Kármán constant, c = ω/k the wave phase speed, U r and θ r the reference wind velocity and direction input at the height z r . The wind velocity U r is the relative wind speed velocity, corrected by the current. It is adapted from Janssen, 1991 with a correction for the friction velocity as in Chen and Belcher, 2000 given by Eq. (2.12) to include the sheltering effect of long waves on short waves, with s u a tuning sheltering coefficient ranging from 0 to 1: (2.11) (2.12) The wave age Z is given by Eq. (2.13) with z α a wave age tuning parameter: (2.13) The roughness height z 1 is evaluated from the following set of equations (2.14): (2.14) The friction velocity u * is evaluated with a law of the wall from the reference wind velocity U r known at the height z r . The roughness height then depends on the wave stress τ w which is evaluated by Eq. (2.15) from Janssen, 1991, and τ is the total stress related to the friction velocity. The roughness height is evaluated by taking into account the effect of high-frequency waves with z 0 given by the Charnock relation and α 0 the Charnock coefficient. It is capped with a user defined value z 0,max to avoid unrealistic high values in extreme high wind conditions. In Eq. (2.15), c and c g are the phase speed and wave group velocity modulus, and χ is the normalized vertical component of the wave-induced velocity in the air, found to be related to the stress tensors: Several tunable parameters are appearing in the parametrization of the wind input described above but the β max parameter is the only parameter globally and homogeneously affecting the wind input source term. It is directly controlling the amount of energy put into the model from the wind input with no distinction on the wind speed, wave direction or frequency. It is controlling the overall wave growth rate and significant wave height bias in the model. Different sets of parameters are offered by default in WAVEWATCH III ® . They were found after a meticulous spectral analysis of the different source terms involved in the package (Ardhuin et al., 2010), reaching a balance between the dissipative terms and the growth terms in key parts of the spectrum notably. However slightly different input and forcing data-sets are used here, bearing different biases. Adjusting the parametrisation is therefore a necessity. As done in Stopa, 2018, we decided to focus only on the growth rate β max .

Input fields
The different input fields that we used are described in this section. The spatial and temporal resolutions of each product are highlighted. The way the data is processed when applicable is also briefly mentioned.

Bathymetry: GEBCO
The latest version of GEBCO is used (Tozer et al., 2019). It is a processed product offering a global bathymetry data on a regular 1/4 arc minute resolution grid, that gathers and merges different sources together. Most of the data for Europe comes from the processed product EMODNET (EMODnet Bathymetry Consortium et al., 2018), which is itself a merged product gathering local surveys from each of the partners in the consortium with a final resolution of 1/16 arc minute. The gaps in the bathymetry are filled with ETOPO1 (Amante and Eakins, 2009), which is a global 1 arc minute resolution processed product using satellite altimetry data.
The bathymetry from GEBCO is then interpolated with a classic bi-linear scheme on each computational grid described previously. The most refined grid is the North-East Atlantic grid with a 3 arc minute resolution, which is a down-grade in resolution compared to GEBCO, meaning the bathymetric product we use is well suited for the wave model presented in this paper.

Wind input: ERA5
The wind growth is the main driving mechanism ruling the generation of waves in our case. As such the wind input data used for the model is the most sensitive input data and its quality and resolution are impacting the accuracy of the model more than the bathymetry or the current field. The ERA5 data-set (Hersbach et al., 2020) is a reanalysis offered by ECMWF using a 4D-Var data assimilation scheme on one of their forecast systems, covering the period from 1979 onwards. It provides hourly time series for a variety of atmospheric variables on a regular 25 arc minute horizontal resolution grid with 137 hybrid sigma/pressure (model) levels in the vertical. The horizontal resolution of the wind input is well suited to generate the swell in the coarser global grid covering the Atlantic (30 arc minute resolution). However the two refined grids used for the present application have a higher resolution than the ERA5 grid (respectively 6 and 3 arc minute). Therefore there is a loss of spatial accuracy due to the resolution of the wind input. More precisely, scales lower than 15 arc minute resolution spatially (around 17 km at a 53N latitude) are not resolved by the wind input data, and therefore are not resolved by the wave model either. The interpolation of the forcing data is done on the fly by WAVEWATCH III ® with a quadratic interpolation in time and a bi-linear interpolation in space.
The wind speed is also adjusted with the current velocity provided to the wave model. It was shown in Bye and Wolff, 1999 and more recently in Renault et al., 2016 that it can significantly impact the momentum transfer from the atmosphere to the water column. In Ardhuin et al., 2012, a more thorough analysis was conducted that looked specifically at the impact of relative wind speed. A better agreement was found when correcting the 10 m wind speed with the ocean surface current. It was also mentioned that using only a fraction of surface current could be a better choice although it was not specifically tested. The literature suggests a clear feedback from the ocean current on the atmospheric boundary layer, surface wind stress and consequently the wave action wind input and dissipation. It is however still unclear what are the exact parameters influencing this interaction and as a result numerical studies seem to have made arbitrary and empirical choices in their parametrization. In our case we are using the surface currents.
The ERA5 data is showing a negative bias for high winds (Ardhuin et al., 2011;Rascle and Ardhuin, 2013). This bias can be partly corrected by tuning adequately the parameters for the wind input, but a more straightforward approach is to directly increase the wind velocity for high winds. This is done by applying the correction (2.16) below, with W the wind speed and W t , c t two constants that can be tuned to set the threshold where the correction starts to apply. This correction is available in WAVEWATCH III ® . The values used in the present application have been suggested by Ardhuin's team in Ifremer: W t = 23 m.s −1 , c t = 1.08 . (2.16)

Current: GLORYS12
The current and sea level fields are of second importance in most large-scale wave models. It is the purpose of this paper to evaluate how sensitive the model is towards these fields. The global ocean eddy-resolving reanalysis GLORYS12 delivered by Mercator Ocean as part of the CMEMS (Fernandez and Lellouche, 2021) is used for these two fields. It covers the years from 1993 to 2019 and is forced using atmospheric data from the ERA5-Interim reanalysis (Dee et al., 2011). It is assimilating along track altimetry data and gives daily averaged time series of the main oceanic variables of which only the horizontal velocities and surface elevation are used for the present application. The data is provided on a global regular 5 arc minute resolution grid with 50 vertical levels well refined at the surface to allow for a consistent representation of the different air-sea processes. The horizontal resolution is comparable with the level of refinement used in the computational grids of this model, but given the time average only the long-term components of the currents and surface elevation are consistently retained. As mentioned in Lellouche et al., 2018, the submesoscales are not captured, ruling out most of the variability in storms, and the strongest component observed in the data is the geostrophic component. Tides are also not included in this model. The first layer of the GLORYS12 data is used at 0.5 m depth, for the surface currents. The current and surface elevation fields are then interpolated on the fly by WAVEWATCH III ® with a linear interpolation in both time and space.

Wind input sensitivity
The WAVEWATCH III ® model described previously is run for the years 2016 and 2017. Five values for the β max parameter introduced in Eq. (2.11) are tested, ranging from 1.60 to 1.80. The default value of β max for ECMWF winds is 1.43, and the maximum value tested and mentioned in the manual is 1.55. However, in Stopa, 2018, the growth rate is tuned for different wind input sources used in a global wave model application. Values between 1.10 and 2.05 are found. In order to validate the model, statistics are computed, averaging over the full duration of the simulation but excluding the first month to avoid unnatural errors due to the model spinning up. The validation is conducted by comparing the model output to wave buoy and weather buoy data, altimetry data and then comparing the accuracy of the model with a reference wave hindcast simulation done for Ireland by Gallagher et al., 2014. The following statistics are used for the validation using the generic notations X and Y for the time series with N records. When the output of the model is compared against an observation, the model data is interpolated in time and space to match exactly the observation: The formulas are adjusted in the case of circular variables to take into account their periodicity. Below the time series X and Y are assumed to be in radians between −π and +π. For the standard deviation the formula given below is not in radians: the value ranges from 0 to ∞ with 0 describing a distribution with no variance at all. The formulas are: mean m r (X) = atan (m(sin(X)), m(cos(X))) , .

Validation against stations
We first compare the output of the model with in-situ observations. The locations of the Acoustic Doppler Current Profiler (ADCP) and buoys are shown in figure 2. They consist in one ADCP (Inishmaan), two wave buoys (Amets Berth A and Amets Berth B) and four weather buoys (M2, M3, M4 and M5). The ADCP was deployed near Inishmaan at the beginning of the year 2017 with a 2 Hz sampling rate. The original motivation for the deployment was to record and analyse high frequency time series for the surface elevation (Fedele et al., 2019). The ADCP raw data is first processed into 12 min time-averaged spectrum time series, and then into the desired mean wave parameters. The two wave buoys Amets Berth A and Amets Berth B are part of the Irish Wave Buoys Network. They measure a wide range of mean wave parameters, amongst which the significant wave height, peak wave direction and mean period T m02 , available with a 30 min sampling rate, are retained. The peak period is also available for those two buoys but the data is of lower quality than the mean period. For consistency the mean period is also retained for the ADCP near Inishmaan. The four remaining weather buoys belong to the Irish Weather Buoys Network. Only the significant wave height, mean wave period T m02 and mean wave direction are available from those four with a 30 min sampling rate as well. All the directions are reported as the direction where the waves are coming from, with a clockwise rotation direction and the origin corresponds to waves coming from the North (nautical convention).
The instruments were not all recording constantly during the two years of hindcast. The model output data is truncated to match accordingly the observation. Spectral output is requested at the location of the buoys in the model every 10 min. The mean wave parameters of interest are then computed from the spectral output to match the exact quantity recorded by the station. To this effect the relative spectra directly computed by WAVEWATCH III ® at each station are transformed into absolute spectra using the relation (2.1). In our case and in order to reduce the amount of variables to include in the analysis, it was decided to focus only on the significant wave height, available at all the stations, on the mean wave period and on the peak direction when available. If not then the mean direction is used. The wave systems are rarely composed of crossing swells in this region, but the mean direction or mean period can lose their physical meaning in the case of a combined wind sea and swell system, making it theoretically more consistent to rely on the peak direction and peak period when possible. The problem is well posed by Portilla-Yandún et al., 2015, where a wave partitioning method is proposed before deriving the integrated parameters. A direct application of a crossing swell/wind sea problem is shown for instance by Breivik and Christensen, 2020 when computing the Stokes drift in such a situation. We are first looking at the distribution of the different variables observed for the validation, shown in figure 3 and 4, and how the growth rate is impacting their shape. Globally a good agreement is found between the model and the observations despite some buoys being located in nearshore areas where the model is not expected to perform well. More specifically, the significant wave height distribution agrees well with the observations at all locations, although a small bias is sometimes noticed for the position of the peak or its value. The same global comment goes for the mean period except for the ADCP at Inishmaan and the M2 buoy where a significant bias in the position of the peak is observed (respectively +2.3 s and −0.5 s). As for the wave direction the distribution is also well predicted except for the ADCP at Inishmaan where the model is showing a strong bias of −15 o . The model has also a tendency to narrow the distribution of the direction, predicting less variability in the direction of the waves (Inishmaan, M2, M5).
Looking at the impact of the growth parameter β max , the differences induced on the shape of the probability density functions are slim, making it difficult to state anything from this analysis. The peak or mean direction is not impacted at all by the growth parameter, and the mean period only marginally. It is consistent with the way the β max parameter is impacting the model, acting homogeneously on all directions and frequencies. A small but noticeable impact is found on the significant wave height as a higher growth rate flattens the distribution by increasing the density of higher wave events. The impact is not a linear shift of the distribution as one could expect from the formulation (2.11). This is due to the non-linear terms appearing in the wave action and, although we are only modifying the input of energy in the system, the redistribution and dissipation are indirectly modified. Increasing the growth rate seems to improve the agreement for some stations (Amets Berth A, Amets Berth B, M4 and M5) and deteriorate it for others (Inishmaan and M2). For M3, it is hard to conclude as the peak value improves but at the same time its location deteriorates.
In order to quantify the differences observed in the distributions, Taylor diagrams are shown in figure 5, completed by RMSE versus bias diagrams, which are a good addition to the compact synthesis offered by the Taylor diagrams for which a constant bias between the data is not picked Figure 3: Probability density functions of the significant wave height (left column), mean period (middle column) and peak direction (right column), for the ADCP at Inishmaan and the two wave buoys Amets Berth A and Amets Berth B, and for the five different values of β max tested. There is a good visual agreement for all parameters for the Amets Berth A and Amets Berth B stations. For the Inishmaan ADCP, the mean period and peak direction are badly captured by the model, because nearshore effects are probably too strong. The influence of β max is visually appearing as small, but consistent between the stations.  up. The probability density function plots give a broader view of the whole performance but they compare the distributions without checking that the events in the model and in the observation record occur simultaneously. As such they can over-estimate the agreement. This seems to be the case for the ADCP at Inishmaan and the two weather buoys M5 and M2 where the errors in direction are significant. For the first two, the respective error in standard deviation was too strong to be plotted on the Taylor diagram without ruining the scale. Overall the same comment can be made that the growth rate impacts noticeably the significant wave height, but only marginally the mean period and not at all the peak or mean direction. Therefore only the significant wave height is retained for the analysis. A recurrent trend is observed with the bias consistently increasing with the growth rate for all stations, as the mean value from the model increases with the growth rate. The standard deviation also increases with the growth rate for all stations, meaning a higher growth rate puts more variability in the model as the distribution flattens and more high wave values are enabled. Surprisingly the correlation is consistently not impacted by the growth rate, and it can be clearly seen on the RMSE versus bias diagram that, depending on the stations, increasing the growth rate improves or not the performance. To be more precise, stations with a negative bias get a better agreement with a higher growth rate increasing the mean value and reducing the RMSE, whereas stations with a positive bias see their agreement deteriorated with a higher growth rate.
It transpires that using stations for the validation does not provide consistent insights about the behavior of the model. This is probably due to the coarse resolution and the nearshore locations of all the stations. It is likely that the errors introduced by the grid resolution on the bathymetry, wind input and propagation scheme are too important at this level and out-range the correction induced by the calibration of the wind input formulation.

Validation against altimetry data
We now compare the output of the model with satellite altimeter observations. This is more suited for large scale models as the observation is spatially averaging the area covered by the satellite beam, smoothing out at the same time the variability associated with small scale processes. The altimeter data used is provided by the Ifremer Laboratoire d'Océanographie Spatiale (LOS) and publicly available on their servers. Significant wave height measurements from nine altimeter missions are extracted, corrected and merged together (Queffeulou, 2013). A sub-square of the North-East Atlantic grid used in our model is used for the validation, also limiting the longitude at −9 • , excluding shallow water points where the bathymetry in particular induces large errors. A satellite crosses the area of interest in a few minutes and because of memory space constraints the model output is requested hourly. Therefore two or three time steps appear in the time interpolation during the same sweep of a satellite, which is not an issue but is worth mentioning.
The distributions on figure 6 are similar to what was observed for the different stations, showing very good agreement between the model and the altimetry data overall. The higher the wave height the lower the probability of encountering this event, with the last 1% of the distribution corresponding to events with a significant wave height above 9 m. The peak value is found around 2.4 m and is well-captured by the simulation in general, although it is sensitive to the β max value. The smaller the β max value, the more the distribution is shifted towards lower values and the narrower the distribution is around the peak position. As the growth rate increases, so does the occurrence of higher wave height events, showing a broader distribution with the peak position shifted more towards higher values.
The model consistently underestimates the small events whatever the value for β max . The reason   A clear bias appears for the waves below 2 m. The tail of the plot is used to infer an optimal β max value, fitting the model for the extreme wave events.
why it is so was not investigated thoroughly but it is believed to be correlated with the area chosen to compute the distributions. The initial growth of the distribution is observed around 0.5 m, while the model captures it around 0.8 m. A similar comparison done on the whole Atlantic domain, not shown in this paper, gives a better agreement for those extremely low wave height events, closer to what is observed in Stopa, 2018 for instance. Some improvement could be made by adjusting the initial growth of waves for this area specifically, but it is not known whether the initial growth parametrization or low wind data is at fault here. The observations based on the probability density are consistent with the Taylor diagram and bias plot shown in figure 7. It can also be shown on the Taylor plot that the correlation is only marginally impacted by the growth rate, which can be the signature for a strong difference in bias solely, or for events with a low occurrence like extreme wave heights. The standard deviation is however significantly impacted. A value of β max between 1.75 and 1.80 seems to capture more accurately the spread of the distribution. Lower values for the growth rate are narrowing excessively the distribution. However for a value of 1.75 the model is showing a positive bias of 5% with a normalised RMSE of 12%. An optimal value of β max to minimise the bias would be reached for 1.60 or even 1.55 (extrapolating), and still giving a normalised RMSE of 11%, which corresponds to 0.39 m for a mean value of 3.5 m. From this analysis the growth rate seems to impact the model mostly in two ways, acting both on the bias and on the dispersion of the distribution, and there is no absolute optimal value that would optimise both parameters.
Using the bias to deduce an optimal value for the growth rate is not believed to be a good strategy. The probability density plot on figure 6 and the Q-Q plot on figure 8 both show that in fact the low bias is obtained by balancing a positive bias for low wave heights with a negative bias for higher wave heights. The Q-Q plot also shows more clearly that none of the runs are capturing correctly the small wave events as they overestimate the wave heights of young wave systems. Little difference is observed between the runs at this level. We thus conclude that the growth rate β max does not impact the initial growth rate of the waves. On the other hand the larger the wave height the more sensitive the coefficient becomes, with a global tendency of the model to underestimate the more extreme events. This non-linear impact explains the discrepancy between the bias and the normalised RMSE.
Increasing the growth rate value enables to improve the agreement for extreme events. An optimal value is found for 1.75, which also gives the best agreement for the standard deviation. We decided to use this value in our model. A similar strategy was pursued by Stopa, 2018, who also makes the argument that extreme events are more harmful and more sensitive when it comes to coastal hazards, impacts on the communities and coastal infrastructure dimensioning.

Comparison with existing work
The accuracy of our model is compared with a relevant existing wave hindcast (Gallagher et al., 2014) also focusing on Ireland. This 34-year wave hindcast was performed using the unstructured version of WAVEWATCH III ® (version 4.11), forced by the ERA5-Interim winds. Several impactful differences between the model used in Gallagher et al., 2014 and our model can be noted. An unstructured grid is used in GA2014 allowing to refine the nearshore region down to 250 m when the most refined regular grid in our model has an average resolution of 3 km. On the other hand GA2014 relies on wave spectra boundary conditions, taken from ERA5-Interim, when our model does not require boundary conditions at all. The wave hindcast done in GA2014 encompasses 34 years  when our own analysis is performed over only two years (2016)(2017). Therefore some variations in the wave statistics are to be expected regardless of the accuracy of each model, as proved by the mean observed values during each period shown in table 2. One important aspect that was not studied is the computational resources needed for each model. For information there are 40 times more nodes in GA2014 than computational points in our model. Despite the differences between the two models, a good agreement is achieved overall with both models showing comparable wave parameters -see table 2. For all stations the significant wave height and mean period predicted by our model show a slightly better agreement with the buoys than GA2014, with a smaller RMSE and higher correlation coefficient. The bias in each model is not consistent. For some stations and some variables, GA2014 can overestimate the wave parameter and our model underestimate it, and the other way around for a different station. As mentioned previously this can be explained by the complexity of the bathymetry impacting the propagation of swell systems or by the accuracy of the wind data impacting the local generation of wind sea. None of the buoys are located in complex coastal features. Amets Berth A, Amets Berth B, M3 and M4 are facing the open ocean. M2 and M5 are respectively in the Irish sea and Celtic sea, which are more protected areas but still in deep enough water so that our model is still covering the location. Moreover none of the buoys are located in shallow water where the depth induced dissipation or wave breaking would dominate. These processes are better captured by a coastal model like GA2014. The four weather buoys M2, M3, M4 and M5 are all located far enough from the coast so that the resolution of the unstructured grid used in GA2014 is still coarse, similar to that of our model.
As for altimetry data, with the significant wave height statistics shown in table 3, a better agreement is obtained for GA2014 with lower bias and lower RMSE, for a similar correlation co-   efficient. It is however ambiguous to conclude anything from those statistics as for both models they are computed differently. The sub-region used in GA2014 includes the Atlantic and Celtic Sea, whereas in our case only the Atlantic region is included, meaning that we are not comparing at the exact same locations. The difference in the recorded significant wave height is the proof of that, with 2.69 m for the analysis conducted in GA2014 and 4.17 m in our case. The area covered in GA2014 includes more shallower and protected areas, where obtaining a good agreement is more difficult. GA2014 seems to outperform our own model. In conclusion we are satisfied with the performance of our WAVEWATCH III ® model. It seems to predict accurately the swell systems propagating to Ireland with an accuracy matching that of a coastal model like GA2014. However the discrepancy in peak or mean direction suggests that the nearshore processes are poorly captured, which is not a surprise given the resolution and the complex bathymetry of the Irish shelf showing high localised gradients. It is also reasonable to assume that the locally generated wind sea is badly represented in the model, as highlighted by the negative bias for young and small waves in figure 6, but it does not seem to impact too much the accuracy of the model at the locations of the stations. Figure 9: Probability density functions of the significant wave height, comparing the run with surface currents and the run without current interactions against altimetry data. A zoom on the distribution tail is shown in the upper right. The impact of currents is small, and no specific trend is noticeable.

Impact of the current field
The WAVEWATCH III ® model validated in the previous section is now used to evaluate the impact of the current field on the accuracy of the model. The two-year hindcast is compared against a reference solution without current interactions.

Global results
We start with a global overview comparing the two runs against altimetry data. The probability density functions shown on figure 9 are very similar with only marginal and point-wise differences. No trend appears, even for extremely low or high events. Looking at the statistics in figure 10, the correlation coefficient is not impacted at all, and the standard deviation only marginally with a 1% difference. The differences in terms of normalised bias and RMSE are also very small, around 1% for both the bias and RMSE. The solution without currents gives the worst agreement with a higher bias and RMSE.
Despite the comments given in the previous section about using buoys for such a large scale simulation we show some of those statistics in figure 11. This is motivated by the idea that current induced refraction is a local process but the focusing or dispersion of the waves can still be observed down-wave, as mentioned by Ardhuin et al., 2012. The significant wave height, the mean period and the wave direction overall show only marginal differences both in terms of distribution and Figure 10: Taylor and RMSE versus bias diagrams for the significant wave height against satellite altimetry data, comparing the run with surface currents and the run without current interactions. Differences are of the order of 1% for the standard deviation, bias and RMSE. No change in correlation is observed. The run with current interactions gives a better agreement.
in terms of biases. The impact of currents is of the order of 1%. Despite being almost marginal for some stations, the influence of currents is more noticeable with a clear separation between the run without currents and the run with currents (Amets Berth A, Amets Berth B, M2 and M4 notably). For Amets Berth A, Amets Berth B and M4 their location in an area known for having stronger currents can explain this observation. There is no consistency however in terms of statistics impacted. For instance the mean direction standard deviation at M4 shows a 3% difference, but only a marginal difference in bias or RMSE for the mean direction, and no difference at all in terms of mean period. It is also noted that the spread either in correlation coefficient, standard deviation error or bias between the stations themselves is more significant than the spreading induced by currents.
The impact of currents observed against either altimetry data or in situ buoys remains small, especially if compared to how sensitive the model is towards the wind input, which makes any conclusion hazardous to draw. There is however a tendency for currents to overall improve the accuracy of the model. Looking at the significant wave height only, the normalised standard deviation is also shown to be reduced due to the impact of currents. It slightly transpires on the distributions (figure 9), as currents seem to narrow the distribution around peak value.

Spatial variability
Looking at the statistics does not give a lot of information on the process as to how and where exactly the small observed differences appear. In this section the spatial variation induced by the currents and the scale at which they impact the wave propagation are analysed for the present application. It was shown in Ardhuin et al., 2017 andMarechal andArdhuin, 2020 that currents have an impact at short spatial scales, between 10 km and 100 km. These two studies focus on two strong boundary currents (Agulhas and Gulf stream respectively). Our case features the North Atlantic Current but it is already weak and diffused in this region, as highlighted in figure 12 that shows the two-year averaged current field. Stable geostrophic gyres still appear. They are eddies detaching from the main current mostly due to the interaction with the bathymetry. Figure 11: Taylor and RMSE versus bias diagrams for the significant wave height (left), peak or mean period (middle) and peak or mean direction (right), for the seven stations used in the model, comparing the run with surface currents and the run without current interactions. Currents do not always improve the agreement, and the impact of currents on the statistics is marginal compared with the spread between all the stations.  There is a clear correlation between those mesoscale geostrophic structures and the propagation of waves, both in terms of significant wave height and direction, as shown in figure 13. The total average difference between the run with current interaction and the run without interaction is shown, both for the significant wave height and the peak direction. The impact of current is still small. The currents induce a variation of 10 cm at best for the significant wave height and 5 o for the peak direction. A stronger impact is observed in nearshore regions and narrow regions like the Channel passage and the Irish Sea. However, the resolutions of our model and of the GLORYS12 model are too coarse to correctly capture nearshore processes so they are ignored. The peak direction is also showing some unexpected large differences in the North-West regions (see the darker red and blue areas in the left plot of figure 13). It corresponds to a somehow strong change in peak direction. The location of this strong shift is slightly different for the two runs. The existence of this front or the reasons why it is being affected by currents are not explained. Weaker impacts of the current are observed for both the wave heights and directions on the Irish shelf. This is mostly due to currents being weaker on the shelf, as tides are not resolved in the GLORYS12 product used. The two-year time average used to derive the statistics also probably cancels out some part of the impact as there is still some variability in the swell direction reaching this region. It would require additional investigation but one would still expect the refraction features to be spatially coherent to some extent and propagate down-wave, along with the swell towards the coast.
The one-dimensional power spectral density functions are shown in figure 14. Following what is done in Ardhuin et al., 2017, they are computed by averaging over time the two-dimensional power spectral density obtained with a Fourier transform of the currents and wave field, then averaging over the direction of the wave-vector to reduce it to one direction. The wavenumbers are expressed in cycles per kilometer, later written cpk in this paper.
Compared to Ardhuin et al., 2017 the currents in our case carry less energy overall, which is expected as the North Atlantic Current is weaker than the Gulf Stream. A nearly constant slope is observed for the currents spectral density, which is expected from the Kolmogorov theory. The resolution of the GLORYS12 product does not allow to capture scales smaller than 10 km (k < 0.1 cpk). However, the wave model interpolates the current field on the computational grid, which explains the difference in slope for the significant wave heights at those smaller scales. Figure 14: Power density spectra (PDS) of the currents and significant wave height with and without current interaction. Solid lines are for the significant wave height, in cm 2 /(cycles/km), and dashed lines are for the currents, in (cm 2 /s)/(cycles/km). The two-dimensional spectra are computed in a box contained between −29 • and −11 • longitude and 46 • and 59 • latitude. They are converted in one-dimensional spectra by averaging over the direction of the wave-vector. Currents impact and explain the wave field at short scales. We expect the PDS decay for the waves to match that of the currents but we are limited by the resolution of the model and GLORYS12 product to realize this observation.
At scales of 100 km and more (k > 0.01 cpk) the currents show no impact. The spatial variability is most likely explained by the propagation of swell systems. For smaller scales the currents add more spatial variability in the wave model. In Ardhuin et al., 2017 the impact of currents was already noticeable at a 100 km scale (k = 0.01 cpk), whereas in our case the threshold is around 50 km (k = 0.02 cpk). The run with currents predicts 20 times more energy at this scale making the currents themselves explain 95% of the spatial variability at this scale.
Focusing on the scale around 50 km (k = 0.02 cpk), the effect of currents is to align the wave spectral density slope with the current spectral density slope. In Ardhuin et al., 2017 currents are stronger and those slopes are actually equal.

Following tracks
In the previous section it was emphasized that the impact of currents is only visible at short spatial lengths. Both global and station time-averaged statistics were also shown to give little to no evidence about the benefits of enabling current interactions. We now look at concrete examples to see in more detail how the interaction appears in the model. One particular track of the satellite Jason2 that crosses a stationary stable geostrophic gyre was chosen for the example, at three different times, as shown in figure 15. The area chosen is off the coast, in deep enough water so that nearshore effects do not impact the quality of the altimetry data. The satellite only takes a few minutes to cover the area shown. Compared to the hourly output from the model the satellite data can safely be assumed to be a spatial cross-section observation of the area at a given time and not a time-series. All the model outputs are time-interpolated to match the time of satellite observation, and for the first time series (2017-02-13) the model output was added a 35 cm negative bias to better match the observations. This was done to remove the errors possibly induced by the wind input and propagation of the swell systems, as only the impact of currents is of interest here, which appears as a signal perturbation at short scales.
A noticeable difference of the order of 20-50 cm can be observed between the run without current interaction and the run with current interaction. The agreement is not always better for the run with currents. Jason2 has an accuracy of 2 cm, the record frequency is 1 Hz corresponding to 1 s sampling period averaging an area of roughly 5 km (varying with the significant wave height). The sampling frequency is way higher than the model output frequency of 1 h and as a result additional wave components with lower time-scales are included in the record, bringing additional short scale variability in the data-set. This variability, looking at the data, is estimated to have an amplitude of the order of 20 cm, which is unfortunately comparable to the current induced variability.
A clear correlation between strong current gradients and locations of enhanced or reduced significant wave height is observed in figure 15, especially for the first row (2017-02-13) featuring a homogeneous wave field across the region. In the second and third rows the wave field is more heterogeneous. A swell system enters the area of interest for the second row (2017-02-23) whereas in the third row (2017-03-14) a swell system leaves the area. Although a more thorough analysis removing the trend induced by propagating swell would be more conclusive, it looks like in all cases the effect of currents is only contained within the zone of strong current gradients and no coherent shapes emerge down-wave of the refraction. This goes against what is observed in Ardhuin et al., 2012, where current effects do have an impact down-wave. The continuous presence of current features and gradients prevents the propagation by inducing new refraction continuously, and the length scale at which those structures would be coherent might be too short to be observed with the resolution of our numerical model. It is also worthwhile mentioning that in Ardhuin et al., 2012 the observation is made in a much more coastal area featuring complex bathymetry, whereas in this case the waves are propagating in deep water.
Focusing now on the current induced effect, the first row shows a text-book example of current refraction. Looking at the first three boxes (blue, red and orange) following the satellite track, perpendicular to the wave propagation, the response of the model is coherent with what is expected from current refraction. The mean wave direction is aligned with the current direction. In the blue and orange side boxes the currents is way stronger than in the red center box. As a result the wave energy is being refracted in this same center box giving a higher wave height than in a case without current interaction. It also seems to agree more with the altimetry observation. The same process and good agreement is also observed for the four boxes in the third row. However disagreement appears for the purple box in the first row. The run with current interaction predicts an increase in wave height when the satellite shows a decrease. The same goes for the red box in the second row, but this time currents in the wave model induce a decrease in wave height when the observation agrees well with a model without current interaction. In those two cases the situation is more complex as currents turn. It is not straightforward to apply the theory but it shows than current interaction can also induce errors in the wave propagation. The occurrences of those events were not studied.

Conclusion
In this paper we presented a WAVEWATCH III ® model forced with the ERA5 winds and currents and sea level from GLORYS12. The purpose of this model is to generate accurate and high resolution swell conditions for Ireland. The model is described and consists of three nested grids. The wind input with the ST4 package is optimised for the specific model set-up in this paper. Disregarding nearshore wave and weather buoys, judged unsuitable for the resolution in the model, we used satellite data derived from altimetry measurements. We found that the best growth rate value for the application shown in this paper is β max = 1.75, which is higher than the default value but still within the range of values found in the literature. The strategy that we followed was to give the best agreement possible for the extreme wave events, despite deteriorating the agreement for the most frequent wave systems. This strategy also enables the wave model to better capture the spread of the distribution.
Discrepancies in the distribution of wave heights are still found. The peak of the distribution is around 2.4 m. The wave model overestimates it by 0.3 m. Extremely low wave heights are badly picked up and underestimated by the model. The initial growth of the distribution is observed around 0.5 m, while the model captures it around 0.8 m. This error seems to be correlated with the area chosen to compute the distribution. A global comparison, not shown in this paper, gives a better agreement for those low wave height events. On the contrary, extreme wave height systems are in general underestimated and the model was optimised in order to correct for this negative bias. A manual correction of the high winds was also conducted in order to correct this negative bias.
We used the model to investigate the effects of currents on the wave propagation over a large scale model, focusing on the region of Ireland which contains a portion of the North Atlantic Current generating mesoscale eddies on the approach of the Irish shelf. Comparing the model against altimetry data, it is found that taking into account the current interaction reduces the error by 1%. In agreement with the literature we also found that currents impact the wave field at very small scales and explain the majority of the wave height variability at those scales, up to 95% below 50 km. Those processes occur at the limit of the scales that can be resolved by our model. The computational grid resolution is indeed 3 km, and 6 km for GLORYS12 which provides the current field. A strong correlation between the impact on the wave field and the presence of the mesoscale eddies is also found.
We studied in more detail the effects of currents on the wave field looking at particular snapshots. The model is found to capture accurately the refraction of waves induced by currents, but some unexpected and unexplained behaviors are observed as well where the model does not follow the behavior observed by the satellite. On the same topic the impact of currents on the significant wave height is found to be of the same order as the variability inherently present in altimetry derived observations, which is evaluated around 0.5 m. This is deemed to be an issue for objectively Figure 15: Following tracks altimetry data compared with model output. The same track at three different times is shown. Maps of the current (middle) and of the significant wave height (right) for the run with currents are also shown at the satellite track time. The black vector arrows represent the current and the red ones the mean wave direction. The impact of currents is not always seen to improve the agreement although some text-book situations are highlighted. The signal recorded by the altimeter is also showing a strong variability, which is matching the effects of currents, making it difficult to compare the data. assessing the gain in accuracy provided by modelling the current interaction at this scale. A better treatment of the altimetry signal would be needed to allows for more relevant comparisons.
The last feature that was not mentioned nor studied is the impact of currents through the wind input term as the wind speed is corrected by the surface currents. This is still a trending topic and using a relative wind speed has been shown to improve the accuracy of the model (Ardhuin et al., 2012, Renault et al., 2016, Ardhuin et al., 2017.