Sea-state contributions to sea-level variability in the European Seas

The contribution of sea-state-induced processes to sea-level variability is investigated through ocean-wave coupled simulations. These experiments are performed with a high-resolution configuration of the Geestacht COAstal model SysTem (GCOAST), implemented in the Northeast Atlantic, the North Sea and the Baltic Sea which are considered as connected basins. The GCOAST system accounts for wave-ocean interactions and the ocean circulation relies on the NEMO (Nucleus for European Modelling of the Ocean) ocean model, while ocean-wave simulations are performed using the spectral wave model WAM. The objective is to demonstrate the contribution of wave-induced processes to sea level at different temporal and spatial scales of variability. When comparing the ocean-wave coupled experiment with in situ data, a significant reduction of the errors (up to 40% in the North Sea) is observed, compared with the reference. Spectral analysis shows that the reduction of the errors is mainly due to an improved representation of sea-level variability at temporal scales up to 12 h. Investigating the representation of sea-level extremes in the experiments, significant contributions (>20%) due to wave-induced processes are observed both over continental shelf areas and in the Atlantic, associated with different patterns of variability. Sensitivity experiments to the impact of the different wave-induced processes show a major impact of wave-modified surface stress over the shelf areas in the North Sea and in the Baltic Sea. In the Atlantic, the signature of wave-induced processes is driven by the interaction of wave-modified momentum flux and turbulent mixing, and it shows its impact to the occurrence of mesoscale features of the ocean circulation. Wave-induced energy fluxes also have a role (10%) in the modulation of surge at the shelf break.


Introduction
Sea level is considered as a key indicator of climate variability and change (e.g. Church et al. 2013;Stocker et al. 2013;von Schuckmann et al. 2018;Oppenheimer et al. 2019), as it integrates the response of different components of the Earth's system (Storto et al. 2019a) interacting at different temporal and spatial scales. Over time scales ranging from centuries to millennia, a leading role is played by very slow and continuous processes such as lithospheric and mantle deformation due to the melting of ice sheets inducing the process of glacioisostatic adjustment (GIA), giving a contribution to sealevel trends (e.g. Peltier 2004;Spada et al. 2006). Sea-level changes at decadal and interannual time scales are due to density and water-mass distribution variations in the ocean, driven by wind, atmospheric pressure, and heat and water fluxes and barystatic sea-level changes through water-mass exchange between the land and the ocean. Wind stress and atmospheric pressure produce, through mechanical stress, a displacement of the water mass involving sealevel variations, due to displacement of the water column. Variations of temperature and salinity due to heat and water fluxes tend to modify the density structure of the water column (steric effect; e.g. Mellor and Ezer 1995;Storto et al. 2019a), which in turn changes the height of the water column. At the regional scales, lateral fluxes also contribute to sea-level variability, adding complexity to sea-level dynamics. As explained in Pinardi et al. (2014), when considering a limited area of the world ocean, mean sea-level tendency is characterized by lateral mass transport fluxes, which makes the regional mean sea-level substantially different from the global estimates.
The accurate estimate of sea-level rise is one of the most important scientific issues that climate change poses, with a large impact for human population (e.g. Lichter et al. 2011;Bonaduce et al. 2016;Storto et al. 2019a) as it was recognized as the main driver for changes in sea level extremes (e.g. Menéndez and Woodworth 2010;Feng and Tsimplis 2014) influencing the non-linear interactions between tide, surges and waves in coastal areas (e.g. Arns et al. 2015;Vousdoukas et al. 2017).
Sea-level rise during the twenty-first century is expected to be larger than during the twentieth century, even if greenhouse gas emissions stopped now (Hu and Bates 2018). On the other hand, large uncertainties are still associated to contemporary sea-level estimates (MacIntosh et al. 2017) and to the future projections (Nowicki and Seroussi 2018). Understanding the processes contributing to sea-level variability and extremes is of major importance to reduce the sources of uncertainty due to the non-linear interactions among the different sea-level components (e.g. Muis et al. 2019).
Recently, Melet et al. (2018) underlined that the effect of waves on the sea-level rise along the coasts of the World's ocean was probably underestimated. In this sense, Dodet et al. (2019) emphasized the strong influence of wave-induced processes (WIPs; e.g Breivik et al. 2015 andStaneva et al. 2017) on coastal sea-level and claimed that this topic deserves higher attention. Arns et al. (2017) underscored the need of considering the changing non-linear interactions between tides, waves, and surges caused by sea-level rise for coastal risk assessments and coastal protection design. Ponte et al. (2019), reviewing the state-of-the-art of coastal sea-level monitoring and prediction, outlined the importance of sea-level observations and the need of accounting also for WIPs within the priorities for the development of optimal and integrated coastal sea-level observing systems.
Ocean waves contribute to modulate the momentum fluxes between the atmosphere and the ocean and waveinduced processes have direct effects on the ocean circulation which can be noticed from the ocean surface to the mixed layer depth, and indirectly even deeper .
Ocean waves can affect water-level changes (e.g. Staneva et al. 2017) through changes to the ocean surface stress, mixing and circulation (e.g. Wu et al. 2019), and WIPs have a major contribution during sea-level extremes (e.g. Mastenbroek et al. 1993;Dietrich et al. 2011;Staneva et al. 2015).
The continuous improvement of the representation physical processes considered by state-of-the-science ocean general circulation models (OGCMs; e.g. Madec 2016) allows nowadays to obtain reliable information about sealevel variability at the different spatial scales considering both long-term trends (e.g. Storto et al. 2019a) and extreme events (e.g. Staneva et al. 2016Staneva et al. , 2017. The Geestacht COAstal model SysTem (GCOAST) is an integrated modelling system developed at the Helmholtz-Zentrum Geesthacht (HZG), which accounts for wavesocean interactions (coupling), and it shows a good skill comparing with in situ and remote-sensing observations (e.g. Staneva et al. 2016Staneva et al. , 2017Wahle et al. 2017;Wiese et al. 2018;Staneva et al. 2018;Schulz-Stellenfleth and Staneva 2019). Different configurations of the GCOAST system were used in previous studies to assess the effect of ocean-wave coupling on the ocean circulation in the North Sea and Baltic Sea (Alari et al. 2016;Staneva et al. 2017;Wu et al. 2019).
The objective of the present work was to assess the contribution of wave-induced processes to sea-level variability and non-tidal residuals (hereafter surge) using a two-way coupled waves-circulation model, based on a high-resolution (∼3.5 km) configuration of the GCOAST system, which considers the Northeast Atlantic, North Sea and Baltic Sea as connected basins. A 9-year period, spanning 2010-2018, was selected to investigate the seastate contribution to sea-level interannual variability and extremes, focusing on the impact of ocean-wave coupling over the different temporal and spatial scales of variability.
The paper is structured as follows: after this introduction, Section 2 describes the ocean-wave coupled modelling system used in this study, the experimental set-up designed to investigate sea-state contributions to sea-level variability, and the observation datasets used to assess the skill of the OGCM simulations. The synergy with observational records and the effect of WIPs on the sea-level variability and surge, both in the open ocean and over the continental shelf areas, are presented in Sections 3 and 4, respectively. The sensitivity of surge to wave-induced processes is described in Section 5. Section 6 summarizes and concludes. relies on the NEMO (Nucleus for European Modelling of the Ocean) OGCM (NEMO v3.6;Madec 2016) with the enhanced implementation of the wave physics (e.g. Breivik et al. 2015;Staneva et al. 2017). The NEMO setup used within the GCOAST system covers the Baltic Sea, the Danish Straits, the North Sea and has a large extent in the North Atlantic ( Fig. 1) at an eddy-resolving spatial resolution of ∼3.5 km, and using an explicit freesurface formulation (Madec 2016;Staneva et al. 2018). The water column is discretized using 50 hybrid s-z* vertical levels with partial cells to fit the bottom depth shape. The OGCM is forced by momentum, water and heat fluxes interactively computed by bulk formulae using the 1h, 0.25 • horizontal-resolution fifth-generation atmospheric reanalyses from the European Centre for Medium-Range Weather Forecasts (ERA5 ECMWF; Hersbach et al. 2020). Atmospheric pressure and tidal potential (e.g. Egbert and Erofeeva 2002) are included in the model forcings. River run-off is considered as a daily climatology based on river discharge datasets, prepared for the GCOAST using data from the German Federal Maritime and Hydrographic Agency ( Ocean-wave simulations are performed using the spectral wave model WAM. The WAM is a third-generation wave model that solves the wave transport equation explicitly without any presumptions on the shape of the wave spectrum. This model represents the physics of the wave evolution for the full set of degrees of freedom in a 2D wave spectrum. A full description is given by WAMDI-Group (1988), Komen et al. (1994), Günther et al. (1992), Janssen (2004), andECMWF (2019). In this application, the WAM Cycle 4.7 runs in shallow water mode, including bottom-induced wave breaking on a model grid situated between 40 • N to 65 • N and −19 • W to 30 • E, with a spatial distribution of φ × λ = 0.03 • × 0.05 • (∼3.5 km). The 2D wave spectra are calculated for 24 directional bands at 15 • each, starting at 7.5 • , and measured clockwise with respect to true north, and 30 frequencies are logarithmically spaced from 0.042 to 0.66 Hz at intervals of f/f = 0.1. The underlying bathymetry is based on the one-minute global General Bathymetric Chart of the Oceans (GEBCO; http://www.gebco.net) topography. The results of the regional wave model are stored every hour.
The driving forces are the ERA5 wind fields at 10 m above the surface (U10; 1-hourly) produced by a dedicated version of the coupled ocean-wave-atmospheric model system IFS Cycle 41r2 4D-Var of the ECMWF (European Centre for Medium-Range Weather Forecasts). The wind products are made available with a spatial resolution of 0.25 • .
Ocean waves influence the circulation through several processes: turbulence due to breaking waves, momentum transfer from breaking waves to currents in deep and shallow water, wave interaction with planetary and local vorticity, Stokes drift and Langmuir turbulence (e.g. Breivik et al. 2015;Alari et al. 2016).
The NEMO ocean model has been modified to take into account the following wave effects calculated by WAM, as described by Staneva et al. (2017) and Alari et al. (2016): Stokes-Coriolis forcing (Wu et al. 2019), seastate-dependent momentum and energy fluxes and waveinduced mixing. A description of wave-induced forcings and their interaction with the ocean circulation is given in Appendix 1.  Craig et al. 2017) which allows numerical simulations to exchange synchronized information representing different components of the Earth system (e.g. Sterl et al. 2012;Wahle et al. 2017;Varlas et al. 2018). In wave-coupled simulations, surface stress, Stokes drift, significant wave height, mean wave period and wave-induced turbulent energy flux fields are passed over from the wave model WAM to the hydro-dynamical model NEMO. WAM receives sea-surface elevation and ice concentration from NEMO. The exchange of surface currents from NEMO to WAM is possible but not used.

Observation datasets
The sea-level simulations from different experiments were first compared against in situ data at the observation positions, as shown in Fig. 1 (red dots).
In situ observations were obtained from the CMEMS INS TAC, which collects near-real-time data from different ocean monitoring systems for the global ocean and the European seas (Copernicus Marine In Situ Tac Data Management Team 2019). Data are continuously updated and quality controlled, as documented at CMEMS INS TAC, 1 which also describes how to freely access the data.
CMEMS provides also a variety of remote-sensing products (Le Traon et al. 2017). In this study, we use the sealevel anomaly (SLA) from the Tailored Altimetry Products for Assimilation System (TAPAS) dataset (e.g. Pujol et al. 2016;Storto et al. 2019b), which provides data without along-track sub-sampling (Dufau et al. 2016) and makes available the corrections applied to the signals retrieved from the altimetric waveforms (e.g. dynamic atmospheric corrections; Carrère and Lyard 2003;Taburet and SL-TAC Team 2019). Here, we consider along-track data from Jason-2 (J2; black lines in Fig. 1), which has a lifetime (2008-2016) that almost overlaps the time window of our experiments.

Experimental set-up
In this section, we describe the design of the experiment performed to investigate the effect of wave-induced processes on the sea-level variability. First, a reference experiment, hereafter referred to as EXP0, was performed during a 9-year period (2010-2018) with NEMO, neglecting the interaction of wave-induced forcings with the ocean circulation. Three other different ocean-wave coupled experiments were carried out varying the wave-induced processes considered during the numerical integration.
To assess the sea-state contribution, an experiment was performed over the same period of EXP0 and considering an ocean-wave coupled configuration, with all three WIPs activated, hereafter referred to as EXP1.
In order to obtain consistent results and quantify the sea-state contribution to sea-level variability, the same initial conditions were used both in EXP0 and EXP1. The numerical simulations were initiated in 2010 with a 'cold start' of the model (e.g. Bessières et al. 2017), which provides only Temperature (T) and Salinity (S) fields. The initial T and S fields for the Atlantic and the North Sea were obtained from the CMEMS FOAM AMM7 model outputs used for BDY forcing (e.g. O'Dea et al. 2012). In the Danish Straits and Baltic Sea, the data were derived from the CMEMS Baltic Sea ocean reanalysis dataset (Hordoir et al. 2015;Pemberton et al. 2017). Both datasets are tri-linearly interpolated on the GCOAST model grid.
To investigate the sensitivity to the different waveinduced processes, two additional experiments were carried during a 1-year period. A first sensitivity experiment, hereafter referred to as EXP2, was performed considering the combined effect of Stokes-Coriolis forcings and waveinduced momentum flux. Following a similar approach, only Stokes-Coriolis forcings and wave-induced energy fluxes were considered in a second sensitivity experiment, hereafter referred to as EXP3.
Studies in the literature underlined the relevance of ocean-wave interaction during storms (e.g. Mastenbroek et al. 1993;Dietrich et al. 2011). Staneva et al. (2017), investigating the effect of wave coupling in the North Sea during 2013, compared wave-coupled and waveuncoupled model integrations with in situ measurement and found significant improvements in terms of surge temporal evolution during the Xaver storm (4th-5th December 2013) when wave-induced processes were considered. In this study, EXP2 and EXP3 were performed over the same year to investigate the sea-state contribution to sea level during a time window characterized by extreme conditions, qualitatively compare the results with those obtained by relevant studies in the literature, and analyze the signature of sea-state contribution over the GCOAST spatial domain.
The experimental set-up used in this study is detailed in Table 1.

Synergy with in situ and remote-sensing observations
In this section, we describe the results obtained in the experiments in terms of the representation of sea-level variability given by the observations. The rows show the names of the relevant experiments, whereas the columns, from 1 to 4, detail the time window, atmospheric forcings (ATM), boundary conditions (BDC) and wave-induced processes (WIPs) considered. Note that in column 4, ALL stands for the combined effect of WIPs: Stoke-Coriolis forcings (STOC), sea-state-dependent momentum (TAUOV) and energy (PHIOC) fluxes Following the approach proposed by Staneva et al. (2017), in situ measurements and numerical simulations were compared considering surge, obtained filtering the ocean tide (through a tidal harmonic analysis; Pawlowicz et al. 2002) from each data source.
A comparison between in situ observations and numerical simulations in terms of total water-level signals is shown in Appendix 2 ( Fig. 13).
A 3-year period, from 2015 to 2017, was selected to compare with in situ measurements: (i) according to the number of stations available during the analysis period; (ii) to analyze the reliability of the sea-level signals in the experiments after several years of integration.
In general, wave coupling enhanced the system to reproduce more accurately the variability of surge signals given by the observations, as shown in Fig. 2. A coherence analysis (Thomson and Emery 2014) was performed to investigate the reliability of the sea-level signal in the experiments, with respect to the in situ observations (Fig. 3). Spectral coherence (Appendix 2) is typically defined as the correlation between two signals as a function of frequencies (or wavelengths if performed spatially) (Bonaduce et al. 2018;Ubelmann et al. 2015;Ponte and Klein 2013;Klein et al. 2004).
In particular, we here focus on the contribution of seastate to surge at the different time scales of variability.
A spectral window between 5 days and 2 h was selected to consider the surge signals both at the low and high frequencies of variability. The results showed that waveinduced processes increased the coherence of the surge signals in EXP1 over the whole spectral window considered, and a significant phase drift reduction was observed at frequencies between 24 and 48 h, compared with EXP0 (Fig. 3). These results can be qualitatively extended also to the comparison with measurements over different periods during 2010 and 2014 (not shown).
In order to further investigate the sea-state contributions, we looked at the error between the surge signals in the observations and the experiments. Here, we focus on the areas over the GCOAST domain where sea-level signals are strongly influenced by the ocean tide (e.g. Davies 1986;Pohlmann 1996) and where tide-surge interaction plays a major role (e.g. Horsburgh and Wilson 2007;Schrum et al. 2016;Arns et al. 2020): i.e. the Southern North Sea (SNS), Northern North Sea (NNS), and Irish Sea (IRS) (black frames in Fig. 1).
In these areas, the surge signals in EXP0 and EXP1 were both significantly correlated (>0.9) with those retrieved from observed records. Looking at the error in the experiments, Fig. 4 shows the differences between the observed surge signals and those in EXP0 (red lines) and EXP1 (blue lines). Large differences, up to 40 cm, were observed considering the reference experiment and the root mean square error (RMSE) ranged between 4 and 6 cm. Nonetheless, the error in the reference experiment was smaller than those obtained by other model-based studies in the literature (e.g. Vousdoukas et al. 2016;Muis et al. 2020),   Fig. 1). The results for EXP0 (red lines) and EXP1 (blue lines) are shown in the spectral window between 5 days (120 h) and 2 h. Left and right panels: spectral coherence and phase (degrees) in the EXPs with respect to the observations, respectively. X-axis is expressed as periods in hours underling the differences between the experimental set-up (e.g. numerical models, atmospheric forcings, horizontal and temporal resolutions) used by the authors and the one used in this study.
In the following part of this section, the results are also presented in terms of sea-state contribution to the reduction of the error (ER). In order to assess the impact of sea-state contribution with respect to the reference simulation, ER is defined as: where k refers to the kth experiment (with k=1). A value of 50% means that the RMSE in the kth experiment has halved with respect to EXP0. The largest impact was observed considering EXP1 over the SNS with an ER of more than 20% compared with the reference simulation (Table 2). Smaller positive ER were observed also in the NNS (14%) and IRS (10%) due to wave coupling.
It is interesting to notice that comparing the experiments with tide-gauge records as monthly means the results were similar in the experiments ( Table 5 in Appendix 2), in contrast with the results observed considering hourly data and underscoring the relevance of data temporal sampling.
The contribution of wave-induced processes was assessed also by a power spectral comparison. The analysis of spectra in a variance preserving form (Thomson and Emery 2014) is shown in Fig. 5. The power spectra of the error clearly show the differences in the experiments. Here, the reduction of the error at the different frequencies (ER spec ) is defined as the percentage decrease of the error with respect to the reference experiment EXP0 (Table 3). Considering EXP1 (blue lines), the impact of wave-induced processes is noticeable down to 12 h, while the impact is weaker at higher frequencies. This was in agreement with The column values show basic statistics of the compared datasets: correlation coefficient (C), RMSE (expressed as cm) and the error reduction (expressed as a percentage; ER) with respect to the reference experiment Lewis et al. (2019), who investigating the impact of wavecoupling in the North West European Shelf, observed an improved representation of the physical signals at semidiurnal frequency due to the coupling. At frequencies between 1 and 5 days, the largest ER spec were observed in the SNS (∼40%). At these temporal scales, significant contributions were noticed also in the NNS (>25%) and in the IRS (>35%). At higher frequencies, between 12 and 24 h, seastate contribution was smaller and ER spec ranged between 10 and 18%, according to the sub-region considered. At the high frequencies (<12 h), the impact of WIPs was weak in the SNS (∼8%) and negligible in the other areas. At the low frequency (1-5 days), the coherence between the surge signals was high (>0.8) in all the experiments. Considering the coherence values at relevant temporal scales, it is possible to compare the results in each experiment, as shown in Table 3.
In particular, a different coherence value observed at the same temporal scale in the experiments provides evidence about the increased (decreased) level of reliability obtained considering (neglecting) WIPs. Considering wave-induced forcings increased the coherence between the surge signals, compared with the reference simulation, and differences between EXP0 and EXP1 can be seen down to temporal scales between 12 and 2 h. At temporal scales between 12 and 36 h, the coherence increased by 20% in the IRS and SNS when considering WIPs, compared with EXP0, while smaller values were observed in the NNS. At high frequencies (<12 h), the coherence was lower than 0.5 (not significant), and results were compared only qualitatively. We argue that also at these temporal scales the coherence in EXP1 was higher than in EXP0 in all the sub-regions considered. At longer time scales (up to 90 days), the coherence between the surge signals was fairly high (>0.8) in both experiments (not shown). This was in agreement with the findings of Piecuch et al. (2019) who comparing observation-and model-based sea-level signals in Northern Europe observed similar coherence amplitudes at the intra-annual frequencies due to the contribution of wind forcings. At these frequencies, wave-induced forcings made the numerical experiment more reliable to the observations, compared with the reference, underscoring their contribution also at the longer time scales. A significant error reduction (30-40%) was observed at frequencies between 30 and 90 days in all the sub-regions considered due to wave coupling ( Fig. 14 and Table 6 in Appendix 2). This was also in line with what was observed by Piecuch et al. (2019), as ocean waves integrate the effect of winds over the fetch across which they blow (e.g. Bonaduce et al. 2019) and contribute to the ocean circulation through wave-induced processes.
In situ records (e.g. tide-gauge) provide valuable information about sea-level variability along the coasts The results for EXP0 (red lines) and EXP1 (blue lines) are shown in the spectral window between 5 days (120 h) and 2 h. Left panel: power spectra of the surge error with respect to the observed data; the spectra are shown in a variance preserving form (cm 2 ). Right panel: Spectral coherence in the EXPs with respect to the observations. X-axis is expressed as periods in hours (e.g. Church et al. 2013), while satellite altimetry data are fundamental for the understanding of the ocean mesoscale dynamics (e.g. Le Traon et al. 2015). In order to assess the synergy with observations also in the open ocean, the results obtained in the different experiments were compared also with remote-sensing along-track data from the Jason-2 satellite mission, over the period from January 2010 to October 2016. In the comparison, we considered remote-sensing data corrected for all instrumental, range and geophysical corrections, except for the Dynamic Atmospheric Correction (DAC; Carrère and Lyard 2003) and the ocean tide (Carrere et al. 2016), before removing the mean sea surface (e.g. Pujol et al. 2016). See also Dinardo et al. (2018) for further details. Figure 6 shows a comparison between Jason-2 and the experiments, at the satellite along-track positions (black tracks in Fig. 1), in terms of annual mean and percentiles. Annual mean values obtained from satellite data had a range of 20 cm, over the period from 2010 and 2016. Annual 95th percentiles were in the order to 1.3 and 1.5 m over the same Columns 2-4: error reduction (ER spec ) with respect to the reference simulation (EXP0) at different temporal scales; values are expressed as a percentage (%). Columns 5-7: spectral coherence (C spec ; 0.8-0.4) in the reference experiment (EXP0). Columns 8-10: as in columns 5-7, but considering the wave-coupled experiment (EXP1); the values show the temporal scale (expressed as hours) at which the coherence falls below 0.8, 0.6, 0.4 period. The sea-level signals in the experiments showed a good synergy with remote-sensing data, both comparing with annual averages and percentiles. It is worth noting that comparing with remote-sensing data the differences between the wave-coupled and reference experiments were not significant. The latter applies also to the comparison over the different sub-regions considered assessing results with in situ data, as shown in Appendix 2 (Fig. 15). This could be related to the repeat cycle of Jason-2 (∼10 days), which could be sub-optimal to observe the sea-state contributions, compared with the hourly sampling of in situ measurements.

Signature of wave-induced processes
In this section, we describe the results obtained in the reference and wave-coupled experiments, looking at the differences in the representation of sea-level extremes. Figure 7 shows the spatial patterns of 95th percentile surge considering EXP0 over the period from 2010 to 2018. In the reference simulation, the surge extremes ranged between 10 and 50 cm, and 8 and 75 cm during JAS (left panel) and OND (right panel), respectively. The largest surge values can be observed in the German Bight, in both seasons, and in the Baltic Sea during OND.
Comparing the previous results with those obtained in EXP1, it was possible to observe the signature of waveinduced processes in the surge signals. Figure 8 shows the 95th percentile surge difference (left panels) between EXP1 and EXP0. Here, it is interesting to notice that sea-state-dependent processes showed their contribution to sea-level variability over the shelf areas, as well as in the open ocean. This was in agreement with the findings of Lewis et al. (2019) who, analyzing the sea-surface height variability over a 3-month period (DJF) during 2017, found similar spatial patterns associated with the minimum and maximum instantaneous differences between two assimilative experiments due to ocean-wave coupling. In the shelf areas (shallower than 200 m), wave coupling enhanced the model simulation and large-scale positive surge differences (>20 cm) were observed in the North Sea (German Bight) and the Baltic Sea during OND, while smaller values were observed in JAS over the same areas. In  Table 4. The spatial variability of the 99.9th percentile surges (not shown) in EXP1 was similar to the patterns observed considering the 95th percentile surge, showing the largest amplitudes in the German Bight of the order of those observed by Dangendorf et al. (2014) in the Southern North Sea over a  centennial period, both during winter (>2.5 m) and summer (>1.5 m). In Fig. 8, the right panels show the relative difference between the wave-coupled and reference experiments, defined as: where EXP k is the surge signal in the the kth experiment (with k=1,...,3). Positive (negative) values of RD mean that the surge signals in the wave-coupled experiment are larger (smaller) than those in EXP0. In this sense, it is possible to notice that over the continental shelf RD is larger than 20%, regardless of the season considered. In the open ocean (e.g. North Atlantic Drift and Bay of Biscay), where the Rossby radius is significantly larger than in the shelf areas of the North Sea (Hallberg 2013), waveinduced processes showed their contribution at spatial scales associated with the mesoscale dynamics. Both negative and positive differences, up to the order of >20%, were observed in the occurrence of cyclonic and anticyclonic eddies. These results were similar to those observed also over other seasons (e.g. DJFM), both in terms of spatial patterns and range of surge differences between the reference and wave-coupled experiments ( Fig. 16 in Appendix 3). A cluster analysis (Forgy 1965;Hartigan 1975), performed considering monthly 95th percentiles in EXP0 and EXP1, underlined that the differences between the experiments observed over the continental shelf areas and in the Atlantic belong to different patterns of variability (not shown).
The different patterns observed over the shelves and in the Atlantic can be explained by looking at the leading contribution of wave-induced processes over the different areas in the ocean and to their non-linear interactions. In the next section, we investigate the sensitivity of surge to different wave-induced forcings to provide evidence about the wave-ocean interaction which mainly contributes to the surge variability.

Sensitivity to wave-induced processes
In this section, we focus on the sensitivity of surge to different wave-induced forcings over a 1-year period (2013).
First, we show the results obtained when EXP0 and EXP1 were considered during 2013. Over this period, the spatial variability of 95th percentile surge in EXP0 (not shown) was comparable with the one observed during 2010-2018, but the range was larger both in OND (5-85 cm) and in JAS (4-75 cm). Larger values were observed looking at 99th percentile surge, up to >1.5 m during OND (not shown).
Looking at the results obtained when an ocean-wave coupled configuration was considered, Fig. 9 shows the surge differences between EXP1 and EXP0 as in Fig. 8, but considering model outputs during the year 2013.
As observed over the whole period considered in this study, WIPs show their largest contribution over the continental shelf areas in the North Sea and the Baltic Sea, and they have a signature also in the Atlantic. In the German Bight, the differences between the two experiments were up to an order of 25 cm. This was in line with the findings of Staneva et al. (2017) who, investigating the effect of wave-induced forcing on the ocean circulation in the North Sea, found similar patterns over the same area during the Xaver storm. In particular, they observed differences of surge peaks of 50 cm between ocean-wave coupled and uncoupled experiments during the storm, in agreement with the results obtained in this study when 99th percentile surge differences were considered (not shown). This last point makes it possible to argue that the spatial variability of the differences between EXP0 and EXP1 reflects the fingerprint of WIPs during the storm that hit the North Sea in December 2013.
The sensitivity to the sea-state-dependent energy flux was investigated in EXP2, neglecting the wave-induced process (Fig. 10). The combined effect of sea-statedependent momentum flux and Stokes-Coriolis forcings has the largest contribution to surge over the continental shelf areas, and the spatial variability of the differences with respect to EXP0 was qualitatively similar to those obtained considering EXP1. An interesting feature was observed moving from the English Channel towards the continental slope in the Bay of Biscay. In this area, a 10% RD was observed during OND, while negative values were obtained in EXP1 during the same period, underlying how sea-statedependent energy fluxes contribute to modulate the surge in the area, due to a modified turbulent mixing of the water column. Similarly, in the Bay of Biscay, larger positive RD values were obtained both during OND and JAS, associated with mesoscale features of the ocean circulation. In the North Atlantic Drift, larger positive RD were observed during OND, compared with those obtained in EXP1, while the opposite was during JAS over the same area. In this area, the results obtained in EXP1 (Fig. 9) and EXP2 (Fig. 10), reflect the spatial scales of the eddies which characterize the in-homogeneity of this branch of the North Atlantic Current (e.g. Garçon et al. 2001).
Looking at the sensitivity to sea-state-dependent momentum flux, Fig. 11 shows the results obtained in EXP3. The modified surface stress plays a key role in the continental shelf areas, while the combination of Stokes-Coriolis forcings and wave-induced energy flux was negligible over these areas. In the North Sea, small RD values (close to zero or negative) were observed here, compared with those in EXP0. At the Danish Straits and in the Baltic Sea, RD values were slightly positive but not comparable with large positive surge values observed in EXP1. Interesting patterns can be noticed in the Atlantic and in the Bay of Biscay, where the combined effect of Stokes-Coriolis forcings and waveinduced mixing clearly show its contribution, associated with the mesoscale variability of the ocean circulation.
A coherence analysis was performed also here to investigate the contribution of WIPs in the wave-coupled experiments at the different spatial scales, with respect to the reference experiment (Fig. 12). Spectral coherence was computed considering the surge signals in the EXPs over an ocean area representative of the North Sea (0 • E, 8 • E; 54 • N, 58 • N). A wavelength window between 7 and 300 km was selected in order to clearly represent the energy content of the surge signal in the sub-domain, given the spatial resolution of the OGCM used to perform the experiments.
The contribution of wave-induced forcings can be noticed both at the large (>200 km) and small scales (>10 km) of variability. Here a low (high) coherence value means that the effect of WIPs considered in the different experiments is large (small) and that their contribution at the different spatial scales is significant (negligible). In EXP1, spectral coherence decreases below 0.8 at spatial scales between 200 and 100 km and it drops to values smaller  (Thomson and Emery 2014) than 0.5 considering wavelengths between 50 and 100 km. The results in EXP2 were comparable with those obtained in EXP1 over the whole spectral window considered, underlying the limited contribution of wave-induced energy fluxes over this area (during 2013) and in agreement with the findings of Staneva et al. (2017). Conversely, spectral coherence in EXP3 was always higher than in EXP1 and EXP2, showing that wave-modified surface stress has a large impact and it shows its contribution both at large and small scales of variability.

Summary and conclusions
Sea-state contribution to sea-level variability was investigated by means of ocean-wave coupled simulations over a 9-year period (2010-2018). Four different experiments were designed by varying the wave-induced processes considered in the numerical simulations. A reference experiment (EXP0) was performed neglecting wave-induced processes computed by WAM. In order to assess the sea-state contribution to sea-level variability, a wave-coupled experiment (EXP1) was performed over the same period of EXP0 considering Stokes-Coriolis forcings, wave-induced momentum, and energy fluxes. The sensitivity of the system to the WIPs was assessed over a 1-year period (2013) considering Stokes-Coriolis forcings combined alternately with wave-induced momentum and energy fluxes in EXP2 and EXP3, respectively.
The experiments were first compared with sea-level observations to assess the reliability of the results in the numerical simulation.
Comparing with in situ measurements, sea-state contributions in the wave-coupled simulation (EXP1) enhanced the system to reproduce more accurately the observation and significant ER (>20%) were observed in the areas of the GCOAST domain where tide-surge interaction has a major influence on the sea-level variability (e.g. the North Sea and the Irish Sea), compared with the reference experiment. In this sense, it is worth noting that the reference experiment showed a good skill compared with other model-based studies (e.g. Vousdoukas et al. 2016).
When assessing the surge signals in the EXPs using power spectra comparison and coherence analysis with the in situ records, the results showed that WIPs significantly contribute (10-40%) to resolve the surge variability, mainly due to a better representation of processes that act at temporal scales up to 12 h. This is in agreement with Lewis et al. (2019), who observed an improved representation of the physical signals at a semidiurnal frequency in North West European Shelf due to ocean-wave coupling.
Significant contributions (up to 40%) were observed also at timescales between 30 and 90 days. This is in line with the results obtained by Piecuch et al. (2019), who found large coherence amplitudes in Northern Europe at the intraannual frequencies between sea level and wind forcings, which in turn affect the sea state.
A high synergy between remote-sensing data and the EXPs was observed over the period covered by the Jason-2 satellite mission (2010)(2011)(2012)(2013)(2014)(2015)(2016), both in terms of annual extremes and mean conditions, even though the differences between the ocean-wave coupled and reference experiments were not significant, probably due to the sub-optimal sampling of satellite measurements.
The spatial signature of WIPs in the surge signals was investigated comparing 95th percentiles in the ocean-wave coupled simulations with those obtained in the reference experiments. Over a 9-year period, the largest surge values in 566 EXP0 were observed in the German Bight (> 70 cm) and the Baltic Sea during OND. Considering WIPs in EXP1 had a large contribution on the surge extremes (∼20%), both over the continental shelf and in the open ocean. Sensitivity experiments performed during the year 2013 showed the contribution of WIPs over the different areas.
During this period, sea-state contribution to the surge observed in the German Bight, during OND, was in line with the findings of Staneva et al. (2017) who investigated the footprint the Xaver storm (December 2013) over the same area. Wave-induced momentum flux had a major contribution over the North Sea and Baltic Sea, while wave-induced energy flux showed a small impact over these sub-domains. On the other hand, wave-modified mixing had a significant contribution (>10%) during OND at the shelf break (e.g. Bay of Biscay), showing that sea-statedependent energy flux modulates the amplitude of surge in these areas through modified turbulent mixing.
In the open ocean, the spatial patterns observed in the North Atlantic Drift and the Bay of Biscay, associated with mesoscale features of the ocean circulation, were driven by the interaction of wave-modified surface stress and vertical mixing.
When assessing the surge signals in the ocean-wave coupled experiment using coherence analysis with the EXP0 in the North Sea, the results showed that WIPs significantly contribute to surge starting from spatial scales between 50 and 100 km. Similar results were observed in EXP2, underlying the small impact of wave-induced mixing over the period and in the area considered. Neglecting waveinduced momentum flux (EXP3) made the surge signals more reliable with the reference simulation, compared with the other ocean-wave coupled experiments, underlying that the interaction between Stokes-Coriolis forcings and wavemodified surface stress contributes to surge both at the large and small scales of variability.
This study aimed to quantify the sea-state contribution to sea level, considering the interannual variability and extremes, at the regional scale using a high-resolution ocean-wave coupled numerical model.
In the future, ocean-wave coupled experiments should be performed over a multi-decadal temporal scale to investigate the sea-state contributions to the non-linear sealevel variations (e.g. Ezer and Corlett 2012;Arns et al. 2017Arns et al. , 2020 and trends.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.

Stokes-Coriolis forcing
The Stokes drift is the drift in the wave propagation direction induced by the motion of surface waves (Stokes 1847). Fluid particle trajectories in water waves are not perfectly circular, mainly due to different speeds of wave crests and troughs (e.g. Staneva et al. 2017), which sets up a difference between the average Lagrangian flow velocity of a fluid parcel and the Eulerian flow velocity (the Stokes drift). As in the case of wave-induced currents, the Stokes drift is influenced by the Earth's rotation and it gives an additional contribution to ocean currents, known as the Stokes-Coriolis force (Hasselmann 1970): where v S is the Stokes drift vector, p is the pressure, τ is the surface stress andẑ is the upward unit vector. The effects of Stokes drift on the advection of tracers (e.g. temperature and salinity) and the mass transport are also considered in the current implementation of the ocean-wave coupling between NEMO and WAM models in GCOAST (e.g. Wu et al. 2019).

Sea-state-dependent momentum fluxes
The presence of waves greatly affects the wind stress, in particular during storms (e.g. Staneva et al. 2017). When waves grow, they absorb momentum from the atmosphere, and as a consequence ocean currents feel less stress. On the other hand, waves release momentum to the ocean when they break. Therefore, the ocean-side τ oc is defined as: where τ a is the atmospheric stress, τ in is momentum extracted by waves from the atmosphere as they grow, and τ db is the momentum released by waves (negative) to the ocean as they mature and break.
Ocean-side stress balances the atmospheric stress only when the input of momentum by wind is balanced by the release of momentum through breaking (fully developed sea).
The momentum flux τ in is enhanced trough the variations of sea-surface roughness (z 0 ) as waves grow, which in turn is related to the friction velocity u 2 * = τ a ρ a : where ρ a is the air density, g is the acceleration due to gravity and α CH is known as the Charnock constant (Charnock 1955). Janssen (1989) assumed α CH not as a constant but as sea-state dependent: whereα CH = 0.006 (see also ECMWF 2019 for further details). The sea-state-dependent roughness can be used to define the wave-modified drag coefficient: where k is the von Kármán's constant. The momentum flux going into NEMO from WAM depends on the wave-modified drag coefficient, which changes the air-side stress and on the ocean-side stress, which, as already mentioned, depends on the balance between wave growth and dissipation ).

Sea-state-dependent energy fluxes
As waves break, they can introduce strong turbulence in the water column (e.g. during a storm; Alari et al. 2016). Several studies demonstrated the importance of wave generated and induced turbulence both at the surface and at depth (e.g. Jones and Davies 1998;Davies et al. 2000;Babanin and Chalikov 2012).
In NEMO, the wave-induced turbulent kinetic energy (TKE) flux depends on the wave energy factor α CB (Craig and Banner 1994), considered as a constant value (α CB = 100) representative of an average between young and mature seas, and therefore irrespective to the sea state.
Observations and numerical model-based studies have shown that α CB is not constant and it actually depends on sea state (e.g. Gerbi et al. 2009;Fan and Hwang 2017). In this sense, Alari et al. (2016) and Staneva et al. (2017), using the WAM spectral wave model estimated the momentum flux from breaking waves source term , showed the variability of α CB in the North Sea and Baltic Sea. This approach was used also in this study to account for sea-state-dependent energy fluxes.

Total water-level signal in the Southern North Sea
In this section, we compare the results obtained in the reference experiment and in EXP1 with in situ measurements in the Southern North Sea in terms of total water level. Figure 13 shows the water level observed at  Arns et al. (2015) and amplitudes larger than 3 m were observed. In particular, wave-induced forcing enhanced the numerical simulation and the reliability of the signals in EXP1 with those of the observations was larger, compared with the reference experiment. The contribution of wave-induced forcings was particularly relevant to capture the maxima of water-level amplitudes gathered from the observation, which on the other hand were significantly underestimated in the reference simulation.

Monthly means surge signals
In this section, we show the results obtained comparing the surge signals in the experiments with those in the observations considered as monthly means. Comparing the experiments monthly means to tide-gauge records the results were similar in EXP0 and EXP1 and RMSE ranged between 3 and 4 cm ( Table 5) . This in contrast with what observed considering hourly data, when significant differences between the experiments were observed, and underlines the importance of data sampling frequency to consider the sea-state contribution to sea-level variability, as already noticed considering remote-sensing data.

Spectral coherence
Spectral coherence is typically defined as the correlation between two signals as a function of frequencies (or wavelengths if performed spatially) (Bonaduce et al. 2018;Ubelmann et al. 2015;Ponte and Klein 2013;Klein et al. 2004). The spectral coherence between the surge signals in the observations (OBS) and in the experiments is defined as follows: where Cr s and S represent the cross-spectral density and spectral density, respectively, of the signals and k refers to the kth experiment.

Error power spectra
A spectral analysis was performed to compare the surge signals in the observations with those of the experiments considering temporal scales between 10 and 90 days. Figure 14 shows the power spectra of the error between the experiments and the observations. At these scales, the largest error reduction (ER spec ) due to sea-state contributions was at frequencies between 30 and 90 days, compared with the reference experiment as shown in Table 5. The results for EXP0 (red lines) and EXP1 (blue lines) are shown in the spectral window between 5 and 90 days. X-axis is expressed as periods in days. The panels show the power spectra of the surge error with respect to the observed data; the spectra are shown in a variance preserving form (cm 2 )

Comparison with remote-sensing data over different sub-regions
The comparison between satellite altimetry data was performed also over the sub-regions considered to compare with in situ data. The results are shown as monthly mean and 95th monthly percentile in Fig. 15. The amplitude of the signals were larger than those observed comparing the different data sources over the GCOAST domain as 95th annual percentiles and annual means. The largest discrepancies were observed in the Irish Sea, while a better agreement between remote-sensing and simulated data was observed in the Southern and Norther North Sea. On the other hand, the differences between the reference and wave-coupled experiments in the sub-regions were not significant, in line with what already noticed comparing the results over the GCOAST domain.

Appendix 3. Signature of wave-induced processes during DJFM and JJAS
The results obtained looking at the surge differences between EXP0 and EXP1 were extended to other seasons. The same analysis was performed also during DJFM and JJAS, as shown in Fig. 16, and the results were consistent with those obtained in OND and JAS (Fig. 8), respectively. Similar patterns of variability were observed both over the continental shelf and in the open ocean areas.