A study of the wave effects on the current circulation in Galway Bay, using the numerical model COAWST

A model coupling waves and currents is set up for Galway Bay, Ireland, using COAWST. The impact of the coupling on the accuracy of the model is shown to be marginal on the overall statistics. The model is in good agreement with the available in-situ observations that are used for the validation. However a closer look at the data during storms shows a deterioration of the agreement, which is barely improved by the coupling. A special analysis of the different processes ruling the ocean circulation and wave propagation is conducted during Storm Hector (2018/06/14). It shows a strong wave-induced surge by 7cm in the back of Galway Bay, and a strong response in terms of currents along the coast of Clare of about 30cm/s. Most of the surge is attributed to the impact of the conservative wave-induced forcing terms, but the main process involved is not identified. The response in currents is partly attributed to the conservative wave-induced forcing terms, but mostly to the wave-enhanced surface roughness.


Introduction
Galway Bay and the Connemara region in the west of Ireland are zones of important economical activity, mostly fishery and aqua-farming. The coast is exposed to the strong Atlantic storms, causing regular flooding and damages to cities and marine infrastructures. Forecasting accurately the local ocean circulation and wave propagation is of interest for the community. The Irish Marine Institute (MI) runs several coastal ocean and wave operational models. One of them, the ocean Connemara model, focuses on the area of Galway Bay. Using these models, the MI regularly carries out specific studies for policy makers and public services. These studies include, for example, assistance in search and rescue missions, tracking the spread of sea lice harming salmonid farms or support to native oyster restoration efforts and oyster farming. The quality of the services delivered and the accuracy of the forecasts depend on the capabilities of the numerical models to represent the oceanographic features and phenomena in the region.
Over the last decade several numerical studies featuring various coupled systems have shown that nearshore processes can be highly impacted by the interaction between currents, waves and winds. Bennis et al. (2020) use the 3D ocean model MARS3D (Lazure & Dumas, 2008) coupled with the third generation wave model WAVEWATCH III ® (The WAVEWATCH III R ® Development Group, 2019), and focus on the Alderney Race (France), a region known for its strong tidal currents. Jakovljević et al. (2018) also study the strong interaction between waves and tides in the Alderney Race, with an application on the tidal energy resource. They use the coupling system COAWST (Kumar et al., 2012;Warner et al., 2008), where the ocean model ROMS (Shchepetkin & McWilliams, 2005) interacts with the wave model SWAN (SWAN Team, 2020). Olabarrieta et al. (2011) reproduce a storm event in Willapa Bay also with the coupling system COAWST.
All these studies highlight a strong interaction of waves and currents. Both surf zone and tidal processes are key to correctly account for the hydrodynamics. In particular tidal currents and the sea level variation are shown to affect wave breaking, which modifies in turn the wave heights and turbulence mixing in the water column. Waves are also refracted by the tidal currents, focusing the wave rays in regions of stronger opposing currents. This last feature is also shown in Ardhuin et al. (2012) where a WAVEWATCH III ® model coupled one-way with a MARS3D model is used for the Iroise Sea (France). A strong focusing is observed near Ouessant and Bannec islands as the waves propagate against a tidal current.
Effects are not limited to surf zone processes. Both in Olabarrieta et al. (2011) and Ardhuin et al. (2012), the local generation of waves is shown to be impacted by the currents. At high frequencies the current-induced Doppler shift can be significant, and the advection of waves can be stopped when facing strong adverse currents. Currents also directly act on the wind stress itself generating the wave growth, where a relative wind is used. However it is not clear how to best account for the currents in this process. Renault et al. (2016) use a coupled air-sea model to highlight a complex retro-effect of ocean currents on the wind in the surface air boundary layer. Waves are also at play in the oceanic surface layer impacting the wind stress momentum flux in the water column. Osuna and Wolf (2005) use the coupled system POLCOMS to investigate the interaction between waves and currents in the Irish Sea. In particular the model features a modified surface roughness, increased by the wave-induced stress, which enhances the surface momentum flux from the wind into the water column. During a storm event they show that the surface stress is significantly increased, leading to a sea level increase up to 15 cm. A similar result is found in Wu et al. (2019), where a similar modification of the surface roughness is done. They are using the ocean model NEMO (Madec & NEMO System Team, 2019) coupled with the wave model WAM (The Wamdi Group, 1988).
The present paper focuses on the possible effects of wave-current interaction on the hydrodynamics of Galway Bay. A few published studies are available focusing either on the currents or waves spawning the Galway Bay area. The main results of those studies will be summarized in Section 2, along with the unpublished knowledge of the area acquired by the MI after years of monitoring and modelling in the region. A fully coupled model using COAWST is set up in Section 3 and described for the area. It is an upgrade of the current operational ROMS standalone model run by the MI. In Section 4, the model is validated over the three periods 2017/01/01-2017/03/27, 2017/05/04-2017/07/28 and 2018/05/10-2018/08/03, highlighting the statistical benefit of the coupling. A more detailed study of the processes introduced by the coupling is then conducted, highlighting the driving mechanisms for this application.

Review of the Galway Bay hydrodynamics
In this section we describe the domain that is used in the numerical model, with a focus on Galway Bay. We also offer a short review of past studies conducted in this area and summarize the main hydrodynamics features highlighted by those studies.
The Connemara region contains several aqua-farms. It features a network of small islands and bays making the hydrodynamics complex. This is highlighted in Jackson et al. (2012) where they Figure 1: Summary map of the oceanographic knowledge of the western Irish shelf, approximately drawn from Nolan (2004). The red dots are tidal gauges, the green dot is the SmartBay location including one ADCP and one wave buoy, the blue dots are the Berth ADCPs (A,B and C from North to South), the yellow dot is the ADCP near Inishmaan. evaluate the cross contamination of sea lice between aqua-farms in Kilkieran Bay using the MI ocean model. However no further studies on the hydrodynamics in the Connemara region has been found.

Irish Shelf
The Irish shelf extends far beyond the limit of the domain studied, shown in Figure 1. Overall the water depths are quite shallow, on average around 200 m and even less (around 100 m) for the domain studied. The waters on the shelf do not interact much with the waters outside the shelf. Further to the West and not shown here is the North Atlantic Current.
A considerable work contributing to the understanding of the hydrodynamics of the western Irish coast is presented in Nolan (2004). A series of hydrographic sections are conducted to investigate the density structure of the water on the Irish shelf and several ADCPs are deployed. They show the existence of a front separating the inland fresh and cold waters and the shelf waters, as shown in Figure 1. The tidal stirring is reduced on the shelf. In summer the wave-induced and windinduced vertical mixing is also weaker. The thermal stratification is therefore able to stabilise the water column and create the front. A jet-like flow is noticed along this front at mid-water depths, from South to North. This Irish coastal current is enhanced by the Shannon River plume. This contributes to increase the mass flux going inside Galway Bay and West of the Aran Islands. Both the location of the front and the Irish coastal current are highly seasonal and more pronounced in summer, again benefiting from the thermal stratification stabilising the water column. The residual currents are mostly wind-driven in autumn, winter and spring, but density-driven in summer, with a location and direction agreeing well with those of the front. The thermal stratification in summer traps cold waters at the sea bottom. A cold dome is present on the Irish shelf, adjacent to the Figure 2: Summary map of the oceanographic knowledge of Galway Bay. The red dots are tidal gauges, the green dot is the SmartBay location including one ADCP and one wave buoy, the blue dots are the Berth ADCPs (A,B and C from North to South), the yellow dot is the ADCP near Inishmaan. front, as shown in Figure 1. Such a cold dome is stable because of the generation of a geostrophic gyre mid-waters, the Coriolis force balancing the density-induced pressure gradient. In theory a decrease in sea level should also be induced at the location of the dome.
In terms of waves, the west coast of Ireland is characterised by rather high significant wave heights, on average around 2.5 m and 3 m (Atan et al., 2016;Gallagher et al., 2014;Gallagher et al., 2016). The coast is regularly hit by strong windstorms and hurricanes over the winter, in general moving from South-West to North-East.

Galway Bay
The Galway Bay area itself has motivated more dedicated research over the years. Its sheltered location and proximity with Galway city make it an attractive area for engineering research projects. It features a marine renewable energy test site .
The bay itself can be divided into an Inner Bay with water depths below 30 m, and an Outer Bay with a maximum depth of 70 m (see Figure 2). The bay is approximately 50 km long and 10 km wide. It is protected from the Atlantic Ocean by the Aran islands. The main connections between the bay and the oceanic waters are the North Sound, north of Inishmore, and the South Sound, east of Inisheer. The main source of fresh water in the bay is through the Corrib River discharge located in Galway city. The bay is also highly influenced by tides, ranging from 5 m to 2 m. It is common for part of the nearshore land to be flooded during strong storms or rain events. Galway Bay is a semi-enclosed basin. As such it prevents most of the swell systems to enter the bay. A wave model consisting of nested grids is set up in Rute Bento et al. (2015) focusing on Galway Bay. It is shown that the dominant wave direction at the SmartBay location is from southwest, which corresponds to swell entering the bay through the South Sound. Wind waves are also shown to be of importance in the bay, with a strong occurrence of short period waves. The model also looses accuracy and underestimates higher wave heights. A poor agreement was found for all wave periods, showing the difficulty of accurately capturing both the propagation of swell inside the bay and the generation of wind waves.
A more accurate analysis of the wave characteristics inside the bay is provided by Atan et al. (2016). Their goal is to report how often optimal conditions for wave energy convectors are met. They conducted a long term analysis over 5 years of the significant wave heights and wave periods at the SmartBay location. They find that the majority of wave systems are characterised by significant wave heights below 1 m and periods below 4 s (44% of events), and 81% of events feature significant wave heights below 2 m and periods below 5 s. Those statistics strongly indicate that the bay is often dominated by locally generated wind waves. Some events with wave heights higher than 3 m are still recorded, which are unlikely to be locally generated given the maximum fetch in the bay. Those events appear on the time series as local peaks. This indicates that the strongest swell systems are still able to penetrate inside the bay. In agreement with Rute Bento et al. (2015), they also find that on average most of the waves inside the bay are coming from the southwest.
Focusing on the ocean currents inside the bay, Ren et al. (2017) compare results from a numerical model to ADCP and radar measurements using different wind forcing conditions. The ADCP data shows that surface currents are stronger and more variable than the depth-integrated currents, sometimes twice as strong. This significant surface shear and high variability are the result of the surface wind stress. Using accurate atmospheric forcing data is therefore critical to accurately capture the currents inside the bay. They show that wind-induced currents compete against tidaldriven currents. As such the currents are more sensitive to wind during neap tides. The surface layer is in general mostly wind driven, with the mid water column being mostly driven by tides. Wind is also shown to greatly explain the weaker northward current component. Ren et al. (2020) investigate the synoptic characteristics of the surface ocean currents measured by a fixed radar and numerically modelled, using a self-organizing map analysis. Both analyses yield similar results. Tidal driven eastwards and westwards patterns are shown to be dominant but only account for 30% of the data. 12 patterns are retained, carrying significant variability. A strong rotation is noticeable in the current field. The main processes known to drive the currents include tides and winds, but also the Corrib River discharge adding a flux of fresh and cold water in the bay as mentioned in Nolan (2004). Such a river discharge is expected to yield a westwards flow following the coast, due to the impact of the Coriolis force.

Numerical model
The modelling system COAWST is used here. It couples the ocean numerical code ROMS with the third generation wave model SWAN. It is a widely used coupling system with a variety of modules. It includes the interaction between waves and currents, which is the main focus of this paper. The implementation presented in Kumar et al. (2012) is used, following the vortex formalism derived in McWilliams et al. (2004). We will introduce all the equations solved in both the hydrodynamics and wave models, mostly referring to Kumar et al. (2012) and Uchiyama et al. (2010).

Ocean model: ROMS
The hydrodynamic code ROMS is a three-dimensional, free-surface and terrain-following numerical model. It solves the Reynolds Averaged Navier-Stokes equations with the hydrostatic and Boussi-nesq approximations in a finite difference framework. A split-explicit time stepping algorithm is used and detailed in Shchepetkin and McWilliams (2005).
The wave-induced conservative terms are implemented following the vortex force formalism of McWilliams et al. (2004). A multi-scale expansion is carried out introducing three components for the hydrodynamic variables, the wave component varying as the gravity waves, the long-wave component varying as the infra-gravity waves generated by non-linear wave interaction, and the current component varying as the longest temporal and spatial time scales. The conservative wave-induced forcing terms appear in the current dynamics balance. The implementation done in COAWST by Kumar et al. (2012) and Uchiyama et al. (2010) simplifies the problem by merging the infra-gravity and current components, also ignoring several higher-order terms in the wave-induced expressions.

3D equations
For the notation, ∇ is the horizontal nabla operator in the (x, y) directions (x is the eastward coordinate and y the northward coordinate), z is the vertical coordinate with e z the corresponding unit vector pointing upwards with 0 the reference mean sea level. The mean sea surface elevation is ζ, −h is the sea bottom level and H = h + ζ the total depth. The time is t, u = (u, v) is the horizontal mean Eulerian velocity, w the vertical component, (u , w ) the corresponding turbulent velocities, (u St , w) are the Stokes velocities, f is the Coriolis parameter, ν is the molecular viscosity. The dynamic pressure p corresponds to the total pressure divided by the reference water density ρ 0 , ρ is the water density and g is the acceleration due to gravity. The full model closure is not reported in this paper to save space, only the most relevant part to the present research. For a more thorough description the reader can refer to Uchiyama et al. (2010) and Kumar et al. (2012). The k − turbulent viscosity model with turbulent eddy viscosity K m described in Warner et al. (2005) is used for the turbulence closure.
The 3D wave-averaged equations read as follows (for simplicity the terrain-following vertical and curvilinear transformation steps are not included): The wave-induced non-conservative forcing terms are contained in F w and D is the horizontal mixing term. The wave-induced conservative terms consist of the vortex force (J, K) and the Bernoulli head K. They are evaluated as follows with H s the significant wave height, k p the peak wavenumber vector, k p its magnitude, and σ p the intrinsic wave peak frequency (the wave parameters are provided by the wave model, SWAN): The total mean sea level ζ (called composite sea level in Kumar et al. (2012)) is decomposed into a dynamical component ζ d and a static componentζ. The static component includes the inverse barotropic effect with P atm the atmospheric pressure at the surface, and the wave set-down forced by the waves labelledζ: . (3.7) The Stokes velocities are given as follows, with U St the horizontal depth-averaged Stokes drift: The bottom boundary condition is not impacted by the waves. The surface boundary conditions at z = ζ read as follows, where P is a wave-averaged forcing term (see Kumar et al., 2012 for its expression): (3.9) (3.10) Depth-integrating the conservation of mass and simplifying with the boundary conditions gives a simple expression for the evolution of the total sea surface elevation, with U the depth-averaged horizontal velocity. The sea level is updated in COAWST using this expression: At the sea bottom z = −h, a logarithmic bottom friction model is used to evaluate the momentum bottom stress τ b = (τ xz b , τ yz b ), not shown here: (3.12)

2D equations
The hydrodynamic code ROMS uses a split time-stepping algorithm. The 2D depth-averaged equations read as follows: The terms are written so that they match the budget output given by ROMS. They consist of the acceleration term (accel), the horizontal advection (hadv) which contains the Stokes drift contribution, the Coriolis force (cor), the wave-induced Stokes Coriolis force (fsco), the surface stress (sstr) and bottom stress (bstr), the pressure gradient (prsgrd), the K vortex force (kvrf), the horizontal J vortex force (hjvf), the non-conservative wave-induced acceleration (fwav) and the horizontal viscosity (hvisc).
The pressure gradient can be further decomposed as follows. The pressure solved in ROMS is the composite pressure p + K, so only the surface Bernoulli head has to be evaluated. Note that the inverse barotropic effect has been separated from the wave set-down effect: . (3.14) The budget terms on the right hand side are also given by ROMS. They consist of the Eulerian sea level adjustment (zeta), the quasi-static sea level adjustment (zetw), the quasi-static pressure gradient (zqsp) and the Bernoulli head contribution (zbeh). The vertical integration is still required to obtain the corresponding 2D budget term. It is straightforward for the wave-induced terms which are only evaluated at the surface.

Non-conservative wave-induced acceleration
The wave-induced non-conservative force terms enabled for this application consist of the acceleration induced by bottom streaming B bf , by whitecapping B wcap and by depth-induced wave breaking B b : (3.15) Bottom streaming The interaction of waves with the sea bed dissipates wave energy in the bottom wave boundary layer. The energy dissipated is transferred to the water column and generates bottom streaming, a wave-induced acceleration aligned with the wave propagation. The wave-induced acceleration appears as the wave orbital velocities become slightly in phase, thus generating an additional stress. The bottom streaming formulation from Uchiyama et al. (2010) is used, with a body force proportional to the wave bottom dissipation bf evaluated by the parameterization of Reniers et al. (2004). It decays in the upward direction following the distribution f bf focused at the bottom, with k wd a decay length related to the bottom wave orbital velocity: (3.16) Whitecapping Whitecapping can occur at any water depth. This wave breaking is caused by a wave steepening of the shortest waves, mostly induced by strong winds. The wave model directly computes the dissipation term wcap . The induced acceleration in the water column is parameterized by the following formulation in Kumar et al. (2012). It mimics the parameterization used for the bottom friction, with f b a distribution function focused near the surface: (3.17) Depth-induced breaking In a similar fashion, the dissipation of energy due to depth-induced breaking b is directly computed by the wave model. The associated induced acceleration in the water column is then computed as follows: (3.18)

Wave-enhanced mixing
Wave breaking of any sort enhances the vertical mixing in the water column (Agrawal et al., 1992). This effect is implemented in Kumar et al. (2012) as a surface boundary condition for the turbulent kinetic energy k, following Feddersen and Trowbridge (2005), and the option is enabled in the model.

Wave-enhanced surface roughness
The bulk parametrization of the surface momentum and heat fluxes involves a surface roughness z 0 which accounts for the sea state. It is traditionally evaluated with the Charnock relation, but it has been shown that a constant Charnock coefficient fails to adequately describe several data sets. The formulation from Drennan et al. (2005) is implemented in COAWST, where the surface roughness is evaluated through the wave age c s,p /u * , with c s,p the wave phase speed at the peak frequency: (3.19) A higher wave roughness increases the surface stress and transfer of energy from the air to the ocean. A rough sea is seen to increase the surface roughness which yields stronger wind-induced surface currents.

Wave model: SWAN
The third generation wave model SWAN (v41.20AB) is used to model the wave propagation with the COAWST modelling system. It solves the wave action balance accounting for current effects. The full derivation of the equations is shown in Holthuijsen (2010). It is worthwhile noting that a vertically uniform current is assumed in the derivation. The wave action balance is given by Eq. (3.20) below, with A(x, y, t, σ, θ) the wave action, (x, y) the eastward and northward horizontal coordinates, t the time, σ the intrinsic wave frequency, θ the spectral direction, c g = (c g,x , c g,y ) the intrinsic wave group velocity, U w = (U w,x , U w,y ) the current field passed down to the wave model, c θ and c σ the rates of change of the direction θ and intrinsic wave frequency σ respectively, k = (k x , k y ) the wave vector with k its modulus, h the mean water depth, ζ the sea surface elevation, and S a global source term: (3.20) (3.21) (3.22) (3.23) The relative variables are all derived from the dispersion relation where the absolute wave frequency accounts for the Doppler shift due to currents: (3.24)

Wave source terms
The source terms used in this application are gathered in S and can be divided into the following expression. S nl is the nonlinear wave interaction term using the standard Discrete Interaction Approximation from Hasselmann and Hasselmann (1985). S in is the wind wave growth input term using the formulation from Yan (1987). S ds,w is the wave dissipation due to whitecapping modelled from van der Westhuysen et al. (2007), S ds,b is the bottom friction (Madsen et al., 1988) and S ds,br is the depth-induced wave breaking (Battjes & Janssen, 1978): (3.25)

Effects of the ocean model on the wave propagation
We briefly highlight here the different effects of the coupling for the wave propagation. Most of those effects have already been mentioned in the introduction.
Total sea level In the previous equations describing the wave propagation the total sea level is adjusted by the mean sea level computed in the hydrodynamics model. In nearshore locations, it induces a significant change, and notably modifies the depth-induced dissipative processes. During low tides more friction and more wave breaking is expected.

Current
The dynamical coupling sends a current field to the wave model. The depth-integrated current of Kirby and Chen (1989) is used, which accounts for the vertical shear in the water column. The wave action balance is derived for a vertically uniform current. With a more realistic sheared current those equations can still be used with U w an effective current field, which in reality is a higher order correction for the wave group velocity. Since this computation is done in the ocean model, an approximation is made by using the peak wavenumber k p : The current field passed down to SWAN mostly impacts the advection through Eq. (3.23). It induces additional refraction through Eq. (3.21) and redistributes energy between the frequency bins through Eq. (3.22). The contribution of currents to wave advection mostly impacts young waves characterised by higher frequencies and lower phase speed and group velocity. Currents can block the wave energy from propagating if they are strong enough to overcome the wave group velocity. The current-induced refraction should impact equally all waves, regardless of their age or frequency. However, the effect is spatially more coherent for narrow-banded swell systems.

Doppler shift
The Doppler shift directly induced by Eq. (3.24) can greatly modify the wave frequencies, especially for locally generated wind waves which are characterised with higher frequencies. This is because of the wave-number multiplying the current term U w : young waves usually have high wavenumbers whereas swell systems have shorter wavenumbers.

Summary of the runs tested
In order to highlight the specific influence of the different processes introduced previously several run configurations are compared in this paper. The most relevant ones are gathered in Table 1. Three time periods are run, from 2017/01/01 to 2017/03/27 (labelled 1701), from 2017/05/04 to 2017/07/28 (labelled 1705) and from 2018/05/10 to 2017/08/03 (labelled 1805). The baroclinic time-step for ROMS is set to 15 s with a 1 s barotropic time-step. The time-step for SWAN is set to 3 min and the same time-step is set for the coupling between the models.

Computational grid
The same computational grid is used for both ROMS and SWAN. The specifications are reported in Table 2. The domain is shown in Figure 1, along with the bathymetry and location of the different stations. Compared to other studies focused on nearshore processes (Bennis et al., 2020;Olabarrieta et al., 2011), the resolution here is quite coarse. However, several studies (Benetazzo et al., 2013;Osuna & Wolf, 2005;Wu et al., 2019) used even coarser domain resolutions without issues, focusing on more offshore processes. The dimensions of the grid are a constraint of this problem. Since there is a possibility for operational applications, the wall-clock time must be kept Runs  low. The grid used here also mimics the domain already used by the MI for their standalone ROMS operational model of Galway Bay. The vertical grid contains 20 levels. A terrain-following σ transformation is applied refining the levels at the surface in order to better capture the wave-ocean and air-ocean interaction at the surface. A coastline file from Ordnance Survey Ireland is used to generate the mask (https: //data.gov.ie/dataset/). The grid is technically curvilinear but given the squared boundary it is very close to a regular grid.

Forcing fields
All the input data used in our model is mentioned here for clarity, along with possible processing that was deemed necessary.

Bathymetry
The bathymetry put in the model is mostly relying on the Integrated Mapping for the Sustainable Development of Ireland's Marine Resource (INFOMAR), a 20 year programme to map Ireland's seabed. It includes a great number of boat and LIDAR surveys. The resolution of those surveys ranges from 1 m to 7 m. Following the best practices for processing bathymetry data reviewed in Sikirić et al. (2009), a reduction step is first conducted by averaging the raw bathymetry data around the computational grid points of the domain.
The available surveys do not cover fully the domain. The gaps are filled with the latest version of GEBCO (Tozer et al., 2019). It is a processed product offering a global bathymetry on a regular 1/4 arc minute resolution grid. Different sources are gathered and merged together in GEBCO. Most of the data for Ireland comes for the INFOMAR surveys already mentioned. They use the global 1 arc minute resolution satellite altimetry data ETOPO1 (Amante & Eakins, 2009) to fill the gaps left by INFOMAR. Using this processed product over ETOPO1 is preferred as it includes slightly more sources than only the satellite data.
Additional processing is conducted in order to merge adequately the different sources and prevent instabilities in the model. Given the tidal range observed, a minimal depth of 5 m is then imposed in the domain, to ensure that the wet cells stay wet in all circumstances. Local and global smoothing is then applied to reduce the pressure gradient errors.

Atmospheric forcing
The atmospheric forcing is crucial either for ocean or wave forecasting. The most refined product available for Ireland is the Met Éireann reanalysis MÉRA (Whelan et al., 2018). It covers the period 1981−2018 and continues in real time with a one year delay. The model runs an atmospheric model and assimilates data from different sources. It is forced at the boundary by ERA-Interim.
The resolution of MÉRA is 2.5 km by 2.5 km, using a curvilinear rotated grid aligned with the main wind direction. Hourly output is available for all the necessary variables needed to force the wave model and the ocean model. Those include the wind at 10 m, the surface air pressure, the air temperature and relative humidity at 2 m, the net short-wave solar radiation and net long-wave radiation, the total precipitation.

Lateral boundary conditions
For SWAN, the lateral boundary conditions come from a multi-grid large-scale WAVEWATCH III ® model described and validated in Calvino et al. (2021). Nested grids refine the model on the Irish Shelf. The WAVEWATCH III ® is forced by the ERA5 winds (Hersbach et al., 2020), and no boundary conditions are necessary as the coarser grid encompasses all of the Atlantic Ocean. The full spectrum is provided every 10 min to capture as best as possible the directional variability in the swell systems entering the domain.
The boundary conditions for ROMS consist of the surface elevation and horizontal velocities. They come from a larger scale North East Atlantic ROMS model (NEA_ROMS) run by the MI covering the Irish Shelf (Nagy et al., 2020). It is itself forced by the GLOBAL_ANALYSIS_-FORECAST_PHY_001_024 model (Chune et al., 2021) for the velocities, the OSU TPX08 model (Egbert & Erofeeva, 2002) for the tides, and the hourly ECMWF High Resolution Forecasts (https: //www.ecmwf.int/en/forecasts/datasets/set-i) for the atmospheric forcing. Output is extracted every 10 min at the exact boundary of the Connemara model.  Table 3: List of the stations where in-situ data is available.

River discharges
Three rivers are included in our application: the River Corrib discharging at Galway City and the Rivers Clarin and Kilkieran discharging in the east end of Galway Bay (see Figure 2). Daily climatology data computed from the records held by the Office of Public Works in Ireland is used for the discharge.

In situ data
For the purpose of validating the model different in-situ observations are used. Their location is shown in Figure 2 by the color dots. They consist of two tidal gauges (red dots) installed in Galway Port and Rossaveal Pier, the three ADCPs Berth A, Berth B and Berth C, all deployed in the East part of the bay (blue dots), one ADCP deployed near Inishmaan (yellow dot), one ADCP at the SmartBay location together with a wave-buoy (green dot). The time periods available for each station, along with their exact location and water depth, are gathered in Table 3. All the data is freely available either directly from the MI servers, or upon request. The tide gauges only record the surface elevation variation from a mean level. The datum used for the tidal gauges in Ireland is the Malin Head vertical datum, which is the mean sea level of the tide gauge at Malin Head in County Donegal. In order to avoid any bias issue linked to the reference level taken when comparing model output and observations the mean value of each time series is subtracted, only keeping the surface anomaly. This is relevant for the tidal gauges but also for the ADCPs, where the surface elevation can also be inferred.
The data available for the three Berth ADCPs and SmartBay ADCP is already processed and consists of the sea surface elevation and velocities for the part of the water column in the range of the ADCP. Additional smoothing is achieved by averaging the signal over 30 min intervals for the SmartBay ADCP and 48 min intervals for the Berth ADCPs, keeping the same sampling rate. This is done to reduce the turbulence captured by the instrument, which does not appear in the ocean model. The ADCP was deployed near Inishmaan in early 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). Both wave measurements and current measurements are available at this location. Raw data is available for this ADCP. It is first processed into 12 min timeaveraged spectrum time series, and then into the desired mean wave parameters. The SmartBay wave buoy measures a wide range of mean wave parameters such as the significant wave height, peak wave direction and peak period.
The instruments were not all recording constantly. Model output data is truncated to match accordingly the observation. The few missing data in the time series are filled by linear interpolation. This is necessary for the tidal analysis to ensure a constant sampling rate throughout the whole series. Spectral output for SWAN and three dimensional hydrodynamic variables for ROMS are requested at the location of the stations every 10 min in the model. The mean wave parameters of interest are then computed from the spectral output to match the exact quantity recorded by the station.

Validation of the model
The validation of the model is conducted by comparing the solution from the standalone runs ROMS-0 and SWAN-0 and the coupled model CPL-0 against available stations. The goal is also to observe the expected improvement coming with the coupling.
The statistical definitions are given in Appendix A. All the directions refer to the nautical convention, and are expressed in degrees. It means that they are reported as the direction coming from, with a clockwise rotation and the origin is coming from the North. The time is always reported in Coordinated Universal Time (UTC).

Sea surface elevation
A constant bias in terms of surface elevation is observed in the model. It reduces the sea level by 20−35 cm everywhere in the domain. This a legacy from the boundary conditions used. It was decided not to fix the boundary conditions as it would also break the consistency between the barotropic velocities and the surface elevation at the boundaries. Instead the mean values are removed from both the observation and the model data. As a result, static or quasi-static effects over duration longer than a few months will not be captured, but it doesn't exclude wave effects characterised by short time scales. However it forces the validation in terms of spectral analysis using a tidal analysis to evaluate how well the model propagates the different tidal components.

Tidal analysis
A tidal analysis covering all the stations recording the surface elevation is conducted. The model is forced at the boundary by the MM, MF, Q1, O1, P1, K1, N2, M2, S2, K2 constituents. However, because of the Rayleigh criteria, the P1 and K2 constituents are removed. The Matlab tool described in Pawlowicz et al. (2002) is used. The observation data sets are linearly filled for missing values when needed. For clarity Table 4 first includes the output of the tidal analysis on the observed data. The errors between the models, either standalone or coupled, and the observations are then shown in Table 5 for both the tidal amplitude and Greenwich phase.
The M2 constituent dominates for all the stations, followed by the S2 and N2 constituents. The remaining constituents are of second importance with amplitudes one order of magnitude smaller. It can be noticed that the two tidal analyses done for Galway Port during different time periods give slightly different outputs with a difference of 4 cm for the M2 constituent, and 14 cm for the S2 constituent. This is expected from conducting a tidal analysis over only three months. The two tidal analyses for the Galway Port tidal gauge and the ADCP deployed near Inishmaan during the 1701 period show very similar tidal constituents, both in amplitude and phase. The   the Galway Port and SmartBay tidal analyses during the 1705 period. The three outputs for the Berth ADCPs show excellent consistency, which is coherent since those three ADCPs were deployed around the same area and during the same period of time.
Looking now at Table 5, there is a clear indication that the wave coupling has a marginal impact on the tidal propagation in the model. The errors are similar for both the standalone and the coupled runs. There is also a good agreement between the model and the observations, especially against the three tidal gauge records with a maximum error of 6 cm for the principal M2 component. Slightly higher errors are observed for the five ADCPs, up to 8 cm, but they remain reasonable.

Storm surges
The sea surface elevation in coastal areas is globally explained by the tidal signal, but non-tidal processes can occasionally greatly impact the sea surface. Severe coastal flooding is often caused by a combination of wind, tide and atmospheric pressure. Accurately capturing those events depends greatly on the accuracy and resolution of the atmospheric forcing used.
The surge signal is extracted by computing the detided signal. Figure 3 shows the surge in Galway Port during Storm Doris (2017/02/23), and the surge at the Berth A location during Storm Hector (2018/06/14). For reference, the output from a ROMS standalone run without any atmospheric forcing is included as well in those figures, only capturing the tidal signal. Despite the detidal processing, residual errors are still left in the signals. This is partly due to the refined temporal resolution of the records (6 min) and the relatively short observation duration (around 3 months). It was still decided not to filter or average any further the signal as it would result in damping the surge in the process.    It is interesting to note that in both cases the run with only tidal forcing already reasonably picks up the surge. We believe that the forcing boundary conditions already capture a good portion of the atmospheric impact on the sea level. It also explains why sea surface anomalies are still observed with run ROMS-T despite not modelling anything else than the tidal propagation. The capability of the model to capture accurately surges seems to depend significantly on the capabilities of the parent model used to generate the boundary conditions. For Storm Doris the model fails to capture accurately the strong surge which reaches almost 1 m, whereas for Storm Hector the model overpredicts the surge at the Berth A ADCP location.
In the case of Storm Doris, ROMS-T predicts a 20 cm surge. This is increased by 9 cm (45% increase) when locally taking into account the atmospheric forcing, and increased by another 3 cm (10%) with the action of waves in CPL-0. For Storm Hector, ROMS-T predicts a 49 cm surge, increased by 10 cm (20% increase) for run ROMS-0, and increased by another 3 cm (5%) with CPL-0. Although the overall performance of the model in each case is different, the local atmospheric forcing and wave coupling consistently increase the amplitude of the surge. The atmospheric forcing in the Connemara model has a higher resolution than in the parent model, which explains a better capture of the local low pressure less spatially averaged, giving a stronger inverse barometric effect. The positive impact of currents will be studied below. It still corresponds to a significant increase in sea level, which is of the same order as that observed in Wu et al. (2019) for a different application in the Baltic Sea.

Velocities
The depth-integrated barotropic velocities are compared between the models and the ADCP measurements. A careful manual integration is conducted for the model data, in order to integrate the same portion of the water column captured by the measurements. Indeed, sea bed deployed ADCPs usually cannot measure the first few meters at the bottom. At every time-step the actual portion within the ADCP measurement range is used to integrate the velocities in the model. In practice it contributes to slightly increase the magnitude of the barotropic velocities computed since the smaller near-bed velocities are truncated.

Overall statistics
The statistical agreement is shown in Figure 4. The model clearly fails to capture correctly the velocities at the location of the ADCP deployed near Inishmaan. For the other four ADCPs the agreement is reasonable. Looking at the correlation coefficients, the eastward components are well reproduced by the model, especially for the three Berth ADCPs. On the other hand poor agreement is found for the northward components. The biases are around 0.5 cm s −1 which is overall small compared to the amplitude of the signals around 20−40 cm s −1 . The currents are mostly tide driven and mean values are centered around 0, so the bias is not the best statistics to evaluate the accuracy of the model. Except for the Inishmaan ADCP, the root-mean-square errors are in-between 3−8 cm s −1 for both components. The poor correlation observed in terms of northward velocities does not reflect on the corresponding biases or root-mean-square errors. Indeed it corresponds to the weaker component of the flow. Comparing with the literature, those errors are of the same order as those obtained in Ren et al. (2017), namely around 4.5−6 cm s −1 , where a different numerical model is used for Galway Bay.
The effect of waves on the currents is overall small. Depending on the station and the component picked, the correlation is either improved or deteriorated by less than 0.5%. A consistent negative bias on the eastward component still appears for all the stations, around −0.3 cm s −1 . A similar bias is also noticed for the northward component of the flow at the SmartBay location only.

Time series
Time series of the depth-integrated velocities around Storm Hector (2018/06/14) are shown in Figure 5 for the three ADCPs Berth A, Berth B and Berth C. Instead of the eastward and northward components, the velocities are shown in terms of magnitude and direction. Those two components clearly highlight the strong tidal dominance on the mass flux.
Outside the time window of the storm, the direction of the flow is well captured especially at Berth C. Small discrepancies appear during the storm, especially at the transition during the flood and ebb. The flow rotation in-between tides shows a chaotic behavior for all stations. For instance looking at the Berth B record, on 2018/06/13 04:00:00 the flow turns clockwise from flood to ebb (direction increases), whereas on 2018/06/13 16:00:00 the flow turns anticlockwise from flood to ebb (direction decreases). It is suspected that the wind driven surface currents are responsible for this behaviour. It is not always accurately captured by the models and contributes to deteriorate the statistics especially on the weaker northward component. The velocity magnitude is rather well captured, but the model tends in general to overestimate the current velocities at its peak, by 10 cm s −1 at worst.
Only during the storm window significant differences are observed between the different runs, both in terms of magnitude and direction. Looking at the time series on 2018/06/14 between 00:00:00 and 18:00:00, it is difficult to assess if the model benefits from the coupling. For instance in the case of Berth A, the slack transition is perfectly captured at 04:00:00 by the coupled model CPL-0 in terms of magnitude when ROMS-0 and ROMS-T show a delay. However, CPL-0 predicts

Wave parameters
Two observations are available when it comes to the wave parameters: the ADCP deployed near Inishmaan exposed to the Atlantic swell, and the SmartBay wave buoy inside Galway Bay. Different sea states and agreements are found for the locations. Only the processed mean wave parameters are available for the wave buoy at SmartBay so spectral partitioning is not possible. In order to separate the contribution from the swell and the local wind wave we loop through the peak period time series. If the value is above a specified threshold then we assume that the swell dominates and all the wave parameters are attributed to the swell component. Otherwise they are attributed to wind waves. The threshold value is arbitrarily picked. A value of 6.8 s is found to work best for the SmartBay wave buoy and 8.0 s for the ADCP near Inishmaan.

Inishmaan ADCP
The relevant wave statistics for the ADCP near Inishmaan are shown in Table 6. The total sea statistics are computed without discriminating against the peak period, and as expected the sea state is dominated by swell (98%). The errors for the significant wave height are around 0.4−0.5 m.
It gives a relative error of 15% with a mean value of 3.1 m recorded by the ADCP, which is deemed reasonable. A worse agreement is found for the peak period and a poor agreement is found for the peak direction. Those values consist of higher order moments of the wave spectrum so this is not surprising. The high peak direction errors around 16−17 • can also be explained by the local nearshore bathymetry, which is badly captured in the model and is strongly responsible for wave refraction. A better agreement is found for the swell dominated sea states than for the wind wave dominated sea states. The errors for the peak period must be compared with the mean values in each case, 11.7 s for the swell and 7.4 s for the wind waves, which gives an average relative error of 6.1% for swell and 6.5% for wind waves. The correlation coefficients are always significantly better for the swell dominated sea states than for the wind waves. It shows that the model captures well the propagation of the swell from the boundary conditions.   Overall there is only a marginal impact of the currents on the waves at this particular location. The mean wave parameters stay unmodified for the total sea. Looking into more details, the highly occurring swell dominated sea states are almost not affected by the currents, but the scarcer wind wave states see a small impact from the currents, which increases the agreement, thus reducing the errors for all three wave parameters. Table 7 shows the statistics for the SmartBay wave buoy for the two runs 1701 and 1705. In both cases a good agreement is found for the significant wave height. The low absolute errors must be put in perspective with the low mean wave height value, 0.94 m, for the run 1701 and 0.71 m for the run 1705. It gives relative errors around 23% for both runs. Looking at the overall impact of the coupling, a slight improvement for all the statistics is observed compared to the standalone run, by less than 1%. No significant impact on either spectral component is observed. Both the swell and wind waves seem to be improved by the same marginal amount. Even during the strong wave event highlighted in Figure 6 the time series do not show any strong impact induced by the coupling.

SmartBay
As opposed to what is observed outside of the bay, the wind wave component gives a better agreement between the model and the observations. The correlation coefficient is slightly better for the wind waves as far as the significant wave height is concerned, but way better for the peak period and peak direction. The peak period and peak direction correspond well to the value of each wave partition, and there is a clear indication that the swell is poorly captured inside the bay. This is especially true for the wave direction. Errors are of the order of 30 • with correlation coefficients close to 0. The agreement is significantly better for the wind waves with errors around 20 • and correlation coefficients between 0.7 and 0.9. The swell is characterised by longer wavelengths and is more sensitive to the bathymetry than wind waves. However the quality of the bathymetry inside Galway Bay is excellent, since the raw surveys are conducted with spatial resolutions of a few meters. It is more likely that the propagation of the swell through the narrow entrance points of the bay is poorly captured by the model. The significant wave height still includes both the swell and wind wave components regardless of which one dominates. The small changes in terms of statistics between the two wave components highlight that the sea state is in general always strongly mixed. The partition between swell dominated and wind waves dominated sea states is indeed close to equal, with 65% for the swell and 35% for the wind waves. This is also highlighted by Figure 6, that shows time series of the significant wave height, peak period and peak direction at the SmartBay location for the run 1705 during a strong wave event. Only looking at the measurements the time peak period shows several jumps between low values around 4 s and higher values around 10 s. The low values match well those of the run CPL-0W where the wave model is only forced with wind and only including wind waves, while the higher values match those of the run CPL-0S only including swell. The time series of the runs CPL-0W and CPL-0S for the significant wave height also indicate that the swell and wind wave components both equally contribute to explain the wave heights. It is also worthwhile noticing that the superposition of the two runs CPL-0W and CPL-0S overestimates the result given by CPL-0. This highlights the heavily non-linear behavior of the wave propagation. Figure 7 shows the wave spectra at the SmartBay location from the run CPL-0 at the three marked instants A, B and C, highlighting once again the strong mixed sea state in the bay. A strong swell component appears in all three cases, with a pronounced southwest origin which would correspond to swell entering through the South Sound. An equally strong wind wave component appears, featuring shorter time periods. Finally a weaker secondary swell component still emerges, with an origin more west-southwest which would correspond to swell entering through the North Sound. The wind wave partition is well aligned with the main swell partition for the first case, but for the two following cases it is miss-aligned as observed previously in Figure 6.

Statistical impact of the sea state on the currents
The only station where measurements of the currents and the wave parameters are both available is the SmartBay station, during the run 1705. A storm is isolated for the study and highlighted in Figures 6 and 8, with a maximum significant wave height of 2.3 m inside the bay corresponding to a swell outside the bay reaching 6 m. During this time period the winds are also strong, around 12 m s −1 , and are mostly south, southwest. Winds certainly have an effect on the upper layer of the water column but since they are mostly constant for this period their effect is not studied.
The global statistics shown in Table 7, and the time series shown in Figure 6 indicate that the wave parameters at the SmartBay location are only marginally impacted by the currents, even during a strong storm event. No further analysis is necessary with this aspect of the interaction and instead we focus on the impact of the sea state on the currents.

Impact of the significant wave height
During the peak of the storm, between slices A and B, the currents are badly predicted by the model. The coupling seems to improve the agreement but a large error is still observed. The impact of the coupling during this event is also quite counter-intuitive. Compared to the standalone run ROMS-0 the coupling reduces the mass transport during the tidal flood (current direction of 260 • ) while it Figure 6: Time series of the significant wave height H s , peak period T p and peak direction θ p during a strong wave event at the SmartBay location. The two runs SWAN-0 and CPL-0 are compared against the measurements highlighting the weak impact overall of the coupling as far as the wave propagation is concerned. The runs CPL-0S and CPL-0W highlight the behavior of the swell component and wind wave component at the SmartBay location. They are shown to equally explain the wave heights meaning the sea state inside the bay is strongly mixed. The three slices A, B and C correspond to the time at which the wave spectrum is plotted in Figure 7.  The standalone run ROMS-0 is compared against the three coupled runs CPL-0, CPL-0S and CPL-0W, highlighting the effects of either swell or wind waves on the coupling. Little impact of the coupling is observed outside of the storm window. During the storm the coupling improves the accuracy of the model but it still fails to capture correctly the velocities. Overall there is a surprising good agreement between the three coupled runs. Discrepancies are only observed for the marked cases B and C.  Table 8: Statistics for the depth-integrated current velocity magnitude for the two runs ROMS-0 and CPL-0, comparing a subset of the data only containing calm sea conditions with a significant wave height lower than 1.2 m, and the complementary subset with stronger sea conditions. There is a clear indication that the currents are better predicted during a calm sea, and badly picked up during strong sea events. The coupling slightly improves the accuracy of the model during those stronger wave events.
increases the mass transport during the tidal ebb (current direction of 80 • ). This effect actually goes against the waves, which are mostly propagating eastward. On the other hand, outside the time window of the storm, almost no impact of the coupling is observed.
In order to quantify this behavior the statistics for the SmartBay ADCP are recomputed discriminating against the wave height (Table 8). The subset labelled calm sea only includes data with wave heights lower than 1.2 m, and the subset labelled storm is the complementary set. Including the full set of data gives an acceptable error of 6.3 cm s −1 that is not improved by the coupling. The subset of calm sea gives a better error of 5.5 cm s −1 that is again not improved by the coupling. This is not surprising as during calm seas the impact of waves is minimal, and most of the data correspond to calm sea events. The storm subset gives a higher error of 10.5 cm s −1 for the run ROMS-0. The error is reduced to 9.6 cm s −1 with the coupling. The correlation coefficient is also improved with the coupling. The error is still high but it indicates that the coupling improves the capabilities of the model under strong sea conditions, which is also when the model is struggles the most to accurately capture the circulation inside the bay.

Impact of the sea partition
For most of the time period shown in Figure 8 the two coupled runs CPL-0S and CPL-0W surprisingly show similar results for the currents. Between the two cases A and B the coupled runs CPL-0 and CPL-0W predict similar mean wave parameters (wind wave dominated), as shown by Figure 6. The interaction terms appearing in the coupling mostly involve the mean wave parameters, namely the significant wave height, peak wavenumber and peak wave direction. One would therefore expect the velocities of CPL-0 and CPL-0W to be similar, which is not the case for case B. Similarly for case C, the mean wave parameters of CPL-0 are similar to CPL-0S (swell dominated), but the responses in terms of mass transport of CPL-0 and CPL-0S are still quite different. This highlights that other parameters come into play, and probably a local analysis is not enough to explain the impact of waves on the currents.
Similarly to what was done previously to highlight the impact of the wave heights on the wave effect on currents, the global statistics for the current velocity are recomputed discriminating against the dominating sea partition (Table 9). The wave effects on the currents seem to be slightly better CPL-0 6.7 0.512 Table 9: Statistics for the depth-integrated current magnitude for the two runs ROMS-0 and CPL-0, comparing a subset of the data where wind waves dominate with a peak period lower than 6.8 s, and the complementary subset where the swell dominates. There is a small indication that the model captures better the currents when the swell dominates and that the coupling improves slightly the solution when wind waves dominate, but the differences remain almost marginal.
captured when the swell dominates the sea state, with an error of 5.9 m s −1 against an error of 6.7 m s −1 for when the wind wave dominates. The coupled run CPL-0, compared to the standalone run ROMS-0, only improves statistics when the wind waves dominate, and slightly deteriorates them when the swell dominates. The statistics are much closer than when discriminating against the wave height, which would indicate that using the dominant sea partition as an indicator is not the best option. Because of the strong mixed sea state the dominant sea partition can fluctuate rapidly between swell and wind waves as well, making the statistics less reliable in the end.

Analysis of the coupling during Storm Hector
In the previous section we compared the model results with measurements. We found that there is a good overall agreement for both the wave parameters and the current velocities. Discrepancies appear for the velocities during strong sea states. A mixed sea state is observed inside Galway Bay. Wind waves are accurately predicted but the model fails to capture correctly the swell component. The coupling is unable to improve the capabilities of the model. A small to marginal impact of the interaction is observed on the global statistics. In this section we look more closely at the impact of the wave-induced forcing on the sea level and currents during Storm Hector. There is also a noticeable impact of currents on wave propagation but it is not as striking, and given that we have some doubts on the capabilities of the model to capture accurately the swell partitions inside Galway Bay, we prefer to focus on the wave effects on the currents. Storm Hector passed through the western part of Ireland between 2018/06/13 and 2018/06/15, from south-west to north-east. At the peak of the storm, wave heights up to 3 m were observed inside Galway Bay.

Wave effects on the surge
The tidal analysis conducted in the previous section (Figure 3, top plot) shows a 3 cm increase of the sea level during Storm Hector induced by the coupling, comparing the run CPL-0 to ROMS-0. Whether or not this contributes to an increase in accuracy is difficult to assess. For instance the bottom plot of Figure 3 shows that the coupling increases the error. However, we argued that most of the observed surge signal is contained in the boundary conditions, so the coupling could very well increase the capabilities of the model by capturing accurately the local wave effects.  Table 1. To save space the runs CPL-3A, CPL-3B, CPL-3C and CPL-4 are not reported but no difference is observed compared to CPL-0.

Averaged surge and wave effects
Assuming now that the model captures accurately the effects of the waves we have a closer look at the surge induced during Storm Hector. The small changes in sea level, of the order of 10 cm at best, are negligible compared to the tidal range. We therefore rely on the difference in sea level between different test runs to observe the wind or wave effects, shown in Figure 9. All the plots correspond to a time average over the peak of the storm, between 2018/06/14 02:00:00 and 2018/06/14 06:00:00. The wave-induced surge is quickly losing intensity, which is why we reduce the window for the time-average, focusing on the strong response at the peak. Plot (A) shows CPL-0 minus ROMS-0, which includes a wide range of wave-induced processes. The difference in sea level is the strongest in the back of Galway Bay and Kilkieran Bay with a maximum increased sea level of 6 cm induced by the coupling.
Plot (B) features ROMS-0 minus ROMS-T and highlights the impact of the atmospheric forcing. It is not shown here but with a better tuned configuration of the model the atmospheric response is attributed solely to the wind impact. Looking at the scaling it can induce a surge up to 30 cm. The wave-induced surge is almost four times smaller than the wind-induced surge, but it is not small enough to be considered negligible and justifies the inclusion of wave effects.
Plot (B) shows CPL-5 minus ROMS-0, which corresponds to the impact of the conservative wave-induced terms. It explains most of the impact of waves up to roughly 4.5 cm. Oddly enough the same wave-induced increase in sea level is observed out of the bay, and also attributed to the conservative forcing terms. The vortex force formalism is composed of different forcing terms. The contribution of each of those terms will be studied below in more depth to identify which effect actually dominates.
Plot (D) corresponds to CPL-0 minus CPL-2 and highlights the impact of the wave-enhanced surface roughness. It explains a smaller portion of the surge up to 1.5 cm, which is significant and should not be neglected. The runs CPL-3A, CPL-3B, CPL-3C and CPL-4 are not reported here but no significant difference is observed. The wave breaking and bottom streaming contributions have a negligible impact on the sea surface elevation.

Explaining the wave-induced surge
The wave-induced surge inside Galway is attributed principally to the conservative terms, and to a smaller extent to the increase in surface stress caused by the wave-enhanced surface roughness (Figure 9, plots C and D). The overall impact reaches 7 cm, which is smaller but not negligible compared to the 30 cm induced by the winds (Figure 9, plots A and B), especially in a populated area like Galway, which is regularly suffering from winter floods during strong storms. We believe that the impact of waves can improve the accuracy of the forecast and the quality of the warning.
The wave-enhanced surface roughness is seen to increase the surface stress. It induces a stronger mass transport towards the coast, in the back of the bay. The momentum terms are introduced and commented later, but one can refer to Figure 14. The enhanced wind surface stress and increased mass transport is met by a stronger pressure gradient, which explains well a strong variation in sea level. Figure 10 shows the contribution of the conservative wave-induced terms. Plot (A) features the surface elevation difference between the runs CPL-5 and CPL-5NOVF, where CPL-5NOVF has the same configuration as CPL-5 except that the vortex force (J, K) is removed from Eqs. (3.1,3.13). The divergence of the Stokes drift appearing in Eq. (3.11) is also removed, as well as the advection by the Stokes drift appearing in Eq. (3.13). Those terms are all evaluated with the Stokes drift depth-integrated mass transport in COAWST. It is difficult to separate them any further because they are evaluated all at once using the Lagrangian velocity within a flux form formalism. Removing those terms provides a good agreement for the overall impact shown in Figure 9, plot (C). It will be argued later, when studying the momentum budget, that the direct impact of the vortex force is marginal, which means that only the mass and momentum transports are responsible for the wave-induced surge.
The small discrepancies are attributed to the Stokes-Coriolis force highlighted on plot (D), showing CPL-5 minus CPL-5NFSCO. The Stokes-Coriolis force induces a small reduction in sea surface elevation in the north-west part of the domain, and an increase in the south-east, consistent with the general direction of the Stokes-Coriolis force pointing south-east.
Plot (B) shows CPL-5 minus CPL-5NOSTK, which stands for CPL-5 only without the divergence of the Stokes drift appearing in Eq. (3.11). This plot highlights the specific effect of the mass transport by the Stokes drift on the sea level. It is strong in Kilkieran Bay and on the east side of the Aran Islands, up to 30 cm. This is balanced by the effect of Stokes drift momentum advection in Eq. (3.13), highlighted in plot (C) with CPL-5NOSTK minus CPL-5NOVF. Unfortunately it was not possible to run a configuration where only the momentum advection by the Stokes drift was removed because the run was too unstable. It is understandable given the strong induced changes in the observed velocity by almost 1 m s −1 .
The direct impact of the divergence of the Stokes drift on the sea surface elevation in Eq. is well explained by an accumulation of water mass transported by the Stokes drift. The divergence is strongly negative near the shore where the incoming Stokes drift is not balanced by an outgoing Stokes drift because of the coastline. As a result it causes an increase in sea level, which is directly linked to the strength of the Stokes drift. This explains why the effect is stronger where the Stokes drift is also observed stronger ( Figure 16, plot C). The Stokes drift magnitude is mostly the result of a reduction in the water depth putting more weight in the surface levels where the drift is more consequent, especially in areas facing the open sea where both swell and wind waves can contribute. However, this is well balanced by the advection of momentum by the Stokes drift. Similarly to the water mass, momentum is also directly advected as shown in Eq. (3.13), with an opposite effect on the sea surface elevation which is not well understood. The response in terms of currents is also well balanced, as shown by the plots (E) and (F). It only features long-shore induced currents. The increase in sea level on plot (B) is met with a clockwise long-shore circulation on plot (E), while the decrease on plot (C) is met with an anticlockwise longshore circulation on plot (F). Long-shore circulations are strongly bounded by the Coriolis force, which means that plot (B) features an offshore current response deviated to its right, coherent with the nearshore increase in sea level. On the other hand plot (C) features an onshore current response, coherent with the increased advection caused by the Stokes drift. For plot (B) the increase in sea level is the cause of the current response, while in plot (C) it is the consequence of it. Both almost cancel each other perfectly, with only a small residual in the two bays.
To sum up, the Stokes-Coriolis force is observed to explain marginally the wave-induced response to the conservative forcing terms. A strong balance between the Stokes drift mass transport and Stokes drift momentum transport is observed, especially in Kilkieran Bay and on the west side of the Aran Islands, which explains most of the response. The Stokes drift mass transport appears as the leading mechanism, with the transport of momentum being the control mechanism. This is understandable as the transport of momentum does not directly impact the sea surface elevation. In both cases it is triggered by the coastline stopping the advection, and forcing a long-shore circulation. No fundamental argument is found as to why the increase is stronger in the back of Galway bay and in Kilkieran bay compared to the rest of the domain.

Wave effects on the circulation inside Galway Bay
The hydrodynamics is complex inside Galway Bay. As mentioned in Ren et al. (2017) the currents are strongly driven by the winds as well as the tides. Strong wave events like Storm Hector highlighted in this paper usually come with strong winds as well. A more complete analysis would also include the impact of wind on the circulation but for the sake of conciseness, we decided to focus only on the impact of waves.

Impact on the depth-averaged transport
The specific impact of the different wave-induced effects on the depth-averaged velocities is studied. Figure 11 shows the depth-averaged and time-averaged current field for different run configurations.
On plot (A) CPL-0 minus ROMS-0 shows the overall impact of all the wave-induced terms in the model. A strong response in terms of currents is observed. It is directed against the overall wave propagation and stays constant during the duration of the storm regardless of the tides. During a tidal flood the coupling reduces the inward flow, while during a tidal ebb it increases it, by at most 30 cm s −1 which is significant for the area. The maximal velocities observed in the model are around 80 cm s −1 . This agrees well with the previous observations made during the validation (Figure 8). The wave effects are mostly attributed to the wave-enhanced surface roughness highlighted on plot (B) with CPL-2, and to a smaller extent to the conservative wave-induced terms with CPL-5 shown on plot (C). The sum of those two contributions indeed agrees quite well with the total observed wave effects. A strong wave impact is also observed in the Kilkieran Bay, and even in the nearshore Connemara region although this part of the domain is not shown here. The wave-induced response is also seen to be directed against the wave propagation, and is induced by the conservative waveinduced forcing terms, highlighted by plot (C).

Momentum balance inside Galway Bay
Time series of the dominant momentum terms in Eq. (3.13) are shown on plot (A) in Figure 12. The depth-integrated momentum terms are integrated horizontally in a squared box encompassing Galway Bay to obtain a general idea of which forces drive the overall circulation inside the bay.
Outside the storm the acceleration terms follow closely the pressure gradient. During the storm the surface stress perturbs this balance but both the acceleration and the pressure gradient are still closely correlated. The pressure gradient is mostly dominated by the tidal motion through the gradient of the sea surface elevation. It is then correctly observed in quadrature with the sea level tidal motion. The Coriolis force has a stable and important impact on the circulation. It appears in phase with the tidal motion described by the surface elevation and in quadrature with the pressure gradient and acceleration. This is expected as the Coriolis force is proportional to the current field itself. The amplitude of Coriolis forcing is roughly half the amplitude of the pressure gradient on the  eastward direction, and of similar amplitude on the northward direction. The Stokes-Coriolis force provides the largest direct wave-induced effects, but appears as a second order driving mechanism, which is expected for a wave-induced effect. All the other wave-induced forcing terms, such as the horizontal or vertical vortex force, the wave-breaking induced momentum or the bottom streaming are completely negligible.
The bottom stress is also noticeable. It is always delayed and opposed to the tidal forcing, which is also expected from a friction force. During the storm it seems that the bottom stress in the main eastward direction is enhanced, almost matching the Coriolis term. Outside the storm window the amplitude of the bottom stress is roughly half the amplitude of the Coriolis force.
The intensity of Storm Hector is well illustrated by the contribution of the surface stress. At the peak of the storm on 2018/06/14 the surface stress is the most important forcing term, especially in the eastward direction where its amplitude accounts for more than twice the amplitude of the Coriolis term. It is balanced by an excess in the pressure gradient term in the opposite direction.
The impact of the coupling on the budget balance is shown in Figure 12. Plot (B) features the time series of the difference in budget terms between the runs CPL-0 and ROMS-0. Plot (C) shows these between CPL-5 and ROMS-0 and plot (D) between CPL-0 and CPL-2. When ROMS-0 is involved in the difference, its contribution for the Stokes-Coriolis force is null. Looking at the scaling, the wave effects roughly induce a 10% change in the budget. The impact of the conservative terms highlighted on plot (C) is twice smaller than the impact of the wave-enhanced surface roughness shown on plot (D). The surface stress is clearly seen to capture most of the wave effects on plot (B). It is induced by the wave-enhanced surface roughness and balanced by the pressure gradient as reproduced on plot (D). The Stokes-Coriolis force is the second wave effect, balanced by the Coriolis force as shown on plot (C). The horizontal advection which is enhanced by the Stokes drift is the last significant wave effect appearing in the budget.
Plot (C) is a little redundant but it shows indeed that the Coriolis force balances the Stokes-Coriolis force, and it highlights the contribution of the Stokes drift in the horizontal advection term. During the tidal flood the advection of currents by the Stokes drift causes a reduction of the eastward horizontal advection, and an increase during the tidal ebb. During the two ebb periods covered by the time series the bottom stress in the eastward direction is seen to be enhanced by the waves on plot (B). It is mostly captured by the conservative terms, inducing an increase of the bottom stress during the tidal ebb and a reduction during the tidal flow.
A small semi-diurnal variability is also observed and induced by the wave-enhanced surface roughness, reducing consistently the bottom friction. Those effects are not explained here. The increased bottom stress is also balanced by the pressure gradient.
Two distinct main wave-induced responses are therefore observed. The first one appears from the conservative wave forcing terms and generates a balance between the Coriolis force and the Stokes-Coriolis force, while the second one is induced by the wave-enhanced surface roughness.

Spatial distribution of the momentum balance
The time series in Figure 12 give an overall picture of balance inside the bay, but they do not accurately describe the local balances of the forces involved. In order to smooth out the variability spawning at short temporal scales the budget terms are averaged over the duration of Storm Hector, between 2018/06/13 23:00:00 and 2018/06/14 10:00:00. Figure 13 focuses on the response to the conservative wave-induced terms by showing the dominant forces involved in the balance. They consist in the Coriolis force on plot (B), which seems to Figure 13: Difference in the significant budget terms, between the runs CPL-5 and ROMS-0, averaged over the duration of Storm Hector. Plot (A) features the pressure gradient, plot (B) the Coriolis force and plot (C) the Stokes-Coriolis force. The balance between the Coriolis and Stokes-Coriolis forces observed in Figure 12 is not reproduced for the strong response along the coast of Clare. Instead a balance between the Coriolis force and the pressure gradient appears. Figure 14: Difference in the significant budget terms, between the runs CPL-0 and CPL-2, averaged over the duration of Storm Hector. Plot (A) features the pressure gradient, plot (B) the Coriolis force and plot (C) the surface stress. Although the highlighted processes are different, the balance between the Coriolis force and the pressure gradient is similar to that observed in Figure 13. be balanced by the pressure gradient shown on plot (A), and not the Stokes-Coriolis force shown on plot (C). On the contrary the Stokes-Coriolis force is well distributed spatially inside Galway Bay. This goes slightly against what Figure 12 shows and will be investigated in a later paragraph.
The main wave effect on the current is induced by the wave-enhanced surface roughness. Figure  14 highlights the strongest budget terms that are modified by including this option. A strong balance between the Coriolis force and the pressure gradient is observed again, shown on plots (A) and (B). Although the effect and process highlighted are different, the conservative wave-induced terms and the wave-enhanced surface roughness both induce a similar response both in terms of current and budget balance. A significant difference is the increase in surface stress shown in plot (C). It is even stronger in the back of the bay where the water depth is lower. It is then well balanced by the pressure gradient specifically in the back of the bay. Figure 12, plot (C), shows that the balance between the Stokes-Coriolis and Coriolis forces is not picked by the difference between the runs CPL-0 and CPL-2, and rightfully so. However a strong response in terms of Coriolis force is still observed here.
It is also worthwhile mentioning that the strong local Coriolis force difference is directly related to the strong response in terms of currents and agrees with its southward direction. The figures show the wave-induced effects, but both the full pressure gradient and the full Coriolis force are tidal dependent. During a flood the tidal flow is going inside the bay and the balance is between a Coriolis force pointing east against a pressure gradient pointing west. It is the other way around Figure 15: Time series of the depth-averaged and horizontally averaged currents and Stokes drift over Galway Bay, on the left is the eastward velocity U and on the right the northward velocity V . The wave-induced currents balance well the Stokes drift on average. during the ebb. What those figures show is the wave effect on those forces, which is similar during both flood and ebb. During a flood the Coriolis force is reduced and the pressure gradient enhanced, and vice versa during the ebb.

Current response to the conservative forces
The Stokes-Coriolis force is the only significant wave-induced forcing appearing in the current velocity depth-averaged budget ( Figure 12). It is small overall compared to the other forcing terms, but dominates the wave-induced response and is balanced by a modification in the Coriolis force. Given the formulations of those two forces, it only highlights a direct balance between the depthaveraged Stokes drift and the depth-averaged Eulerian velocity, which is shown in Figure 15, where only the response to the conservative terms with CPL-5 is considered. Overall and on average inside the bay, the Eulerian response to the wave forcing matches quite well the Stokes drift. The amplitude of the response is quite small, around 2 cm s −1 at best in the eastward direction and 1 cm s −1 in the northward direction, but it corresponds to a spatial average. Looking at Figure 13 it is not obvious whether a local equilibrium or at least partial equilibrium is reached between the Coriolis and Stokes-Coriolis forces. The Stokes-Coriolis force is seen quite weak in the center of Galway Bay behind the Aran Islands, and stronger in the inner bay and along the coast of Clare, which also corresponds to the areas where the Coriolis force is stronger. Plot (A) in Figure 16 shows the resultant of those two forces. The areas still left with a strong forcing are well correlated with areas where a response is observed (Figure 11, plot B), and equivalently where the wave-induced Coriolis force is seen strong (Figure 13, on plot B). They are highlighted in the figure by a red ellipse and correspond to Kilkieran Bay, the north side of Inishmore, along the coast of Clare and off to the west of Black Head, and inside Inner Galway Bay. In general those two forces do not locally balance each other, indicating that the Eulerian response is a local feature, which balances the global input of momentum provided by the wave forcing.
The Stokes-Coriolis force appears then as a natural candidate to explain the wave-induced impact on the currents, but it is found unable to capture fully the response as shown on plot (B) in Figure 16. However the full response is well reproduced when adding the impact of the mass and momentum transports by the Stokes drift. Plot (C) shows the sum of the individual responses from the Stokes-Coriolis force with CPL-5NOSFCO and from the vortex force plus the momentum transport induced by the Stokes drift with CPL-5NOVF. Those last two effects are paired due to the implementation done in COAWST, but given the study of the budget terms, the momentum transport is assumed to be the only relevant mechanism and the vortex force negligible. Figure 16: Plot (A) shows the resultant of the Coriolis force plus the Stokes-Coriolis forces, between the runs CPL-5 and ROMS-0, averaged over the duration of Storm Hector. It is well correlated with areas of strong responses observed in Figure 11, plot (B). Plot (B) shows the impact of the Stokes-Coriolis force only with CPL-5 minus CPL-5NOFSCO, while plot (C) shows the impact of the vortex force and Stokes drift mass transport with CPL-5 minus CPL-5NOVF. Both effects are able to individually capture the response to the conservative terms, but added together they explain most of it.

Current response to wave-enhanced surface roughness
Although the response in current to the wave-enhanced surface roughness is similar to the response to the Stokes drift forcing (Figure 11), it doesn't involve the same processes. In terms of balance (Figure 12), the overall forces involved are the surface stress, which is enhanced by the waves, and the pressure gradient balancing the surface stress. Most of this balance occurs inside Inner Galway Bay (Figure 14, plots A and C), and explains well the wave-induced surge. It is highlighted in Figure 17, plot (A).
A strong balance between the Coriolis force and the pressure gradient appears along the coast of Clare on plot (B). However, it is mostly a consequence of the southwards response in current. The Coriolis force response pointing west is induced by the south response in current. It induces a stronger mass transport and a slight water level increase to the west of the current response. A pressure gradient is thus generated that is opposed to the gradient in surface elevation that is emerging. This is verified on plot (C), which shows the difference in sea surface elevation, where the scale has been modified to highlight this particular location. The gradient in surface elevation can be roughly estimated. A variation of 0.7 cm over 2 km gives a pressure gradient of 35 · 10 −6 m s −2 . This shows that the Coriolis force and pressure gradient are able to maintain the balance, but it does not explain the origin of this excess of momentum.
Plot (A) in Figure 18 shows the time series of the depth-averaged eastward and northward velocities averaged inside the green box shown in plot (B). The southward velocities, mostly defining the response along the coast of Clare inside the red ellipse in plot (B), match quite well the eastward velocities mostly strong in the blue box. The strong southward response can hence be explained through the conservation of mass of the enhanced wind-induced velocities due the wave-enhanced surface roughness. The excess of transport going inside Inner Galway Bay has to be compensated by an outward flow, which is contained on the south part of the bay. This is consistent with the global rotation forced by the Coriolis force inside the bay. Plot (C) shows the averaged circulation forced by the wind by showing ROMS-0 minus ROMS-T. It already induces a southward flow along the coast of Clare highlighted by the ellipse. Therefore the wave forcing increases this particular feature of wind-induced circulation. Again, the main circulation is still forced by the tides, so the response just described above consists in a reduction of the inward flood tidal flow, or an increase of the outward ebb tidal flow. In both cases it is induced by a stronger wind-stress caused by the Figure 17: Plot (A) shows the resultant of the pressure gradient plus the surface stress computed for the run CPL-0 minus CPL-2. The contribution of those two budget terms inside the inner bay is well canceled. Plot (B) similarly shows the resultant of the pressure gradient plus the Coriolis force. The contribution of those terms along the coast of Clare is mostly canceled as well. Plot (C) zooms on the sea surface elevation induced by the wave-enhanced roughness, which is estimated to be strong enough to generate the observed response in pressure gradient. The red ellipses are highlighting the areas of interest in each plot. Figure 18: Plot (A) shows the time series of the depth-averaged and spatially averaged eastward and northward velocities inside the green box shown in plot (B). The southward velocities match well the eastward velocities, preserving the conservation of mass. Plot (B) shows the response in terms of currents to the wave-enhanced surface roughness. The eastward transport in the blue box is believed to be balanced by the southward transport in the red ellipse going around Black Head and then through the South Sound. Plot (C) shows the wind-induced circulation with ROMS-0 minus ROMS-T, already featuring a southward response in the red ellipse. action of waves inside Galway Bay.

Conclusion
A coupled wave-ocean model is set up for Galway Bay using COAWST. The model domain is the same as currently used in an operational forecasting system based on ROMS at the Marine Institute, Ireland. The impact of wave-ocean coupling was tested in this study. This coupling comes with a computational and financial cost. It is adamant to assess the impact of the coupling on the capabilities of the model. We follow closely the configuration presented in Kumar et al. (2012). The wave-induced surface streaming is however not included. It makes our model highly unstable, despite several refinements of the surface layers and reduction of the time-step. Candidly tested during Storm Hector, it is observed to have a noticeable impact on the solution, but negligible compared to the wave-enhanced surface roughness for instance.
The model is validated using tidal gauges, ADCPs, a wave buoy, and covering three time periods of roughly four months each. It gives an overall good agreement. The errors in terms of velocities are of the same order, if not slightly less, than a previous model found in the literature for Galway Bay (Ren et al., 2017). Excellent agreement on the tidal signal is found for the surface elevation. The swell outside Galway Bay agrees well with the measurements. It validates the boundary conditions. However a loss of accuracy is observed inside Galway Bay when it comes to the swell partition. This is especially true when looking at the peak wave direction and peak frequency associated to the swell systems. The locally generated wind waves are still very well captured by the model. We suspect that the model fails to capture accurately the depth-induced processes as the swell passes through the North and South Sounds. We also show that the sea state inside Galway Bay is characterised by a mixed sea between one wind wave partition and two swell partitions, entering through the North and South Sounds.
The impact of the coupling is not seen to improve significantly the statistics at the stations where measurements are available. First, the coupling is only significant inside Galway Bay during strong sea conditions. This is necessary to have a strong wave forcing impacting the currents and sea level, but also necessary to have enough wave energy inside the bay that can be significantly impacted by the currents and sea level. The long-term statistics are smoothing out those rare events. The model is significantly less accurate during those storm events, especially for the currents. The coupling is only slightly improving the agreement. An interesting follow-up would be to specifically study the influence of the wind conditions on the accuracy of the model. Strong sea states and storms are commonly accompanied with strong winds that are known to drive significantly the circulation in Galway Bay. Second, focusing on Storm Hector, we show that although a small increase in sea surface elevation is observed everywhere in the domain by 2 cm on average. The strongest effects are observed in the back of Galway Bay, reaching 7 cm. The wave effects are even more local in terms of currents. A weak Eulerian response is observed in the surface layers everywhere in the model, going against the wave propagation (not shown in this paper), but the strongest response in terms of currents appears along the coast of Clare. The response induces a reduction of the tidal flood and an increase of the tidal ebb by 30 cm s −1 at most, which is significant.
The response in sea level is divided between the direct impact of the Stokes drift inducing a strong mass transport directed in the back of the bay, and the wave-enhanced surface roughness increasing the wind-induced momentum and therefore mass transport also pushing water in the back of the bay. The Stokes drift accounts for most of the effect, around 4.5 cm, while the increased wind transport by caused by the waves only accounts for 1.5 cm. In terms of currents the response is also divided between the conservative and non-conservative effects. The Stokes-Coriolis force and the Stokes momentum transport induce a noticeable Eulerian response in the bay. Although the wave forcing is well spread inside the bay, the response is well localised inducing changes in the currents up to 10 cm s −1 along the coast of Clare. When integrated inside the bay it is still found to balance the Stokes drift. The wave-enhanced surface roughness overall increases the wind-induced circulation in the bay. It dominates over the conservative wave effects by accounting up to 20 cm s −1 along the coast of Clare. putational facilities and support. Finally, the authors are grateful to Teledyne RD Instruments (TRDI) for providing the ADCP and for their useful contributions.
Credit authorship contribution statement C Calvino: Conceptualisation of this study, Methodology, Software, Data curation, Writing -Original draft preparation. T Dabrowski and F Dias: Overall supervision -Contributed to manuscript preparation and approved the submitted version.

A Statistical parameters
The following statistical parameters are used for the validation, using the generic notation X and Y for the time series with N records: root-mean-square error : Pearson correlation : 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 values range from 0 to ∞ with 0 describing a distribution with no variance at all: mean : m r (X) = atan (m(sin(X)), m(cos(X))) , standard deviation : σ r (X) = −2 ln m(cos X) 2 + m(sin X) 2 , root-mean-square error : Pearson correlation : . (A.8)