Drivers of Interannual Sea Level Variability on the Northwestern European Shelf

Sea level on the northwestern European shelf (NWES) varies substantially from year to year. Removing explained parts of interannual sea level variability from observations helps to improve estimates of long ‐ term sea level trends. To this end, the contributions of different drivers to interannual sea level variability need to be understood and quanti ﬁ ed. We quanti ﬁ ed these contributions for the entire NWES by performing sensitivity experiments with a high ‐ resolution con ﬁ guration of the Regional Ocean Modeling System (ROMS). The lateral and atmospheric boundary conditions were derived from reanalyses. We compared our model results with satellite altimetry data and used our sensitivity experiments to show that nonlinear feedbacks cause only minor interannual sea level variability on the shelf. This indicates that our experiments can be used to separate the effects of different drivers. We ﬁ nd that wind dominates the variability of annual mean sea level in the southern and eastern North Sea (up to 4.7 ‐ cm standard deviation), whereas the inverse barometer effect dominates elsewhere on the NWES (up to 1.7 ‐ cm standard deviation). In contrast, forcing at the lateral ocean boundaries results in small and coherent variability on the shelf (0.5 ‐ cm standard deviation). Variability driven by buoyancy ﬂ uxes ranges from 0.5 ‐ to 1.3 ‐ cm standard deviation. The results of our sensitivity experiments explain the (anti)correlation


Introduction
Sea level on the northwestern European shelf (NWES) has been rising since the nineteenth century , and this rise is projected to accelerate in the coming centuries (e.g., Church et al., 2013;Oppenheimer et al., 2019). Fluctuations of sea level on various time scales complicate detecting sea level trends from observations, and sea level variability will continue to be the dominant source of uncertainty of sea level change for the coming decades (Palmer et al., 2018). Removing explained parts of variability from sea level records can lead to more accurate estimates of trends and accelerations (e.g., Calafat Thompson, 1986;Wahl et al., 2013). Therefore, the processes driving interannual to multidecadal sea level variability on the NWES need to be understood.
Drivers of sea level variability on the NWES have mainly been studied using the long tide gauge (TG) records that are available in the region. Wahl et al. (2013) found that interannual sea level variability is largest at TGs in the southeastern North Sea. Using linear regression, Dangendorf et al. (2014) showed that a large part of the interannual sea level variability at TGs in the southeastern North Sea can be explained by variability of wind stress, and at the western U.K. and Norwegian coastlines by the inverse barometer (IB) effect (Stammer & Huttemann, 2008). Additionally, a link has been found between the North Atlantic Oscillation (NAO) and sea level on the NWES (Chen et al., 2014;Tsimplis & Shaw, 2008;Tsimplis et al., 2005;Wakelin et al., 2003;Woolf et al., 2003).
On annual to multidecadal time scales, sea level at several locations on the NWES is significantly correlated with sea level along the eastern boundary of the North Atlantic. Miller and Douglas (2007) and Woodworth et al. (2010) correlated sea level at Brest (France) and Cascais (Portugal) to sea level pressure (SLP) over the subtropical North Atlantic and suggested a link between multidecadal sea level variability and mass redistribution related to gyre-scale adjustments. Others showed that sea level variability is coherent from northwestern Africa to the NWES and strongly correlated with integrated longshore winds (LSWs) on time scales of 1 year to decades (Calafat et al., 2012;Sturges & Douglas, 2011). This suggests an important role of boundary dynamics and, in particular, wind-driven coastally trapped waves (CTWs). Dangendorf et al. (2014) hypothesized that CTWs also affect sea level variability in the North Sea. Using the correlation between sea level in the North Sea and the eastern North Atlantic, Frederikse, Riva, Slobbe, et al. (2016) closed the sea level budget for the North Sea for 1958-2014. More recently, Chafik et al. (2019) argued that steric anomalies in the North Atlantic are communicated to the NWES through Ekman transport and affect interannual to decadal sea level variability, since sea level, along-slope winds, and subpolar gyre properties covaried during 1993-2016. These studies have explained sea level variability on the NWES at interannual to multidecadal time scales using statistical relationships between observed sea level at different locations and between sea level and a variety of oceanic and atmospheric drivers. Many of these drivers, such as SLP, winds, gyre circulation, precipitation, and the NAO, are correlated themselves (Woodworth et al., 2009). This makes the correlations between observed sea level and sea level drivers challenging to interpret. Additionally, since TGs are sparsely distributed and located at the coast, TGs alone do not allow us to characterize sea level variability spatially and determine whether observed signals are purely coastal or originate in the open ocean. Satellite altimetry has a better spatial coverage than TGs but can be less accurate in coastal zones due to altimetric corrections and land contamination.
Instead of using observations, Roberts et al. (2016) ran sensitivity experiments with a global ocean model to link interannual sea level variability to different modes of climate variability. More recently, Tinker et al. (2020) used a high-resolution regional ocean model to explore the influence of atmospheric forcing and lateral ocean boundary conditions on interannual sea level variability on the NWES. However, they used a climate control simulation and did not investigate the effect of specific drivers. Hence, the contributions of different drivers to interannual sea level variability in the satellite altimetry era are not yet fully clear for the entire NWES region.
Here, we investigate different drivers of interannual sea level variability for a period in the satellite altimetry era . To this end, we developed a high-resolution setup of the Regional Ocean Modeling System (ROMS; Shchepetkin & McWillams, 2005) for the NWES. We derived boundary conditions for our model from ocean and atmosphere reanalyses (GLORYS12v1 and ERA5, respectively). By separately varying the lateral ocean boundaries, SLP, winds, and buoyancy fluxes in our model, we quantify the effects of these drivers on interannual sea level variability.
We explain our ROMS setup and methods in section 2 and show that our model simulations agree well with satellite altimetry data in section 3. We show the influence of forcing at the lateral ocean boundaries, SLP, winds, and buoyancy fluxes on interannual sea level variability and discuss correlations between variability at different coastal locations in section 4. We end with a discussion in section 5 and our conclusions in section 6.

Methods and Data
Here we introduce our ROMS setup for the NWES (section 2.1) and the reanalysis data that we prescribed at the boundaries of ROMS (section 2.2). Additionally, we discuss the satellite altimetry data that we used to evaluate our model (section 2.3) and present the sensitivity experiments that were used to study the influences of different drivers of interannual sea level variability (section 2.4).

ROMS Setup for the NWES
ROMS is a primitive-equation, curvilinear, hydrostatic, Boussinesq, and free-surface ROMS (Shchepetkin & McWillams, 2005). It runs on a staggered Arakawa C grid and has stretched terrain-following vertical coordinates (Song & Haidvogel, 1994). Terrain-following coordinates are particularly suited for modeling regions with large bottom topography gradients such as the slope from the deep ocean toward the NWES.
Our model domain extends from 36°N to 62°N and from 20°W to 10°E and has four lateral open boundaries ( Figure 1). Our model has 30 vertical levels and a horizontal resolution of (1/8)°by (1/8)°(~11 km), so resolves eddies in the open ocean. With this resolution, the English Channel and Irish Sea are at least four to five grid points wide. The bathymetry (Figure 1), derived from ETOPO1 (Amante & Eakins, 2009), was smoothed with grid stiffness rx 0 = 0.13, which reduces horizontal pressure gradient errors (Beckmann & Haidvogel, 1993). We used mixing along geopotential surfaces.

ROMS Boundary Conditions From Reanalysis Data
We derived surface and lateral boundary conditions for ROMS from reanalysis data for 1993-2018. For the lateral ocean boundaries, we obtained monthly mean values of sea surface height (SSH), zonal and meridional currents, and temperature and salinity from GLORYS12v1 (Fernandez & Lellouche, 2018). GLORYS12v1 is a global ocean reanalysis with a horizontal resolution of (1/12)°by (1/12)°(~8 km) and covers 1993-2018. Its model component NEMO (Madec & NEMO Team, 2016) is driven by atmospheric forcing from the ERA-Interim reanalysis (Dee et al., 2011). GLORYS12v1 assimilates satellite sea surface temperature, in situ temperature and salinity profiles, and along-track satellite altimetry data corrected for the IB effect and the 18.61-year lunar nodal cycle (Fernandez & Lellouche, 2018).
We imposed Chapman conditions for the free surface (Chapman, 1985), Flather radiation conditions for barotropic velocities (Flather, 1976), and mixed-radiation nudging for baroclinic velocities and ocean temperature and salinity (Marchesiello et al., 2001;Orlanski, 1976). We nudged temperature and salinity toward outgoing nudging at a time scale of 1 year and ingoing nudging at a time scale of 1 day, using a 20-point relaxation zone at the lateral boundaries of the domain. Baroclinic velocities were nudged at the boundaries only. To suppress computational noise, we also applied a 20-point sponge layer, multiplying the diffusion and viscosity by a factor that linearly increases to 10 at the lateral boundaries.
ROMS and NEMO (GLORYS12v1) are both Boussinesq ocean models, meaning that local steric effects are modeled prognostically but non-Boussinesq steric effects are not (Griffies & Greatbatch, 2012). Although non-Boussinesq steric effects have a negligible imprint on monthly or longer regional sea level patterns, they do impact global mean sea level (Greatbatch, 1994;Griffies & Greatbatch, 2012). The global mean sea level in GLORYS12v1 was corrected by adding the global mean steric effect to the SSH diagnostically . This signal may enter our model via the lateral boundary conditions but has only small interannual variability.
The atmospheric boundary conditions were derived from ERA5 (Copernicus Climate Change Service (C3S), 2017). ERA5 is based on the IFS41r2 Earth system model and is available for 1979 to a few months before present. The atmospheric forcing was applied using bulk formulae based on (Fairall et al., 1996). We prescribed solar radiation daily using an idealized diurnal cycle and precipitation 12-hourly. Winds at 10 m above sea level, atmospheric pressure, cloud coverage, air temperature, and relative air humidity at 2 m above sea level were prescribed on a 6-hourly basis. We also applied a climatological river run-off based on the river discharge data set of Dai (2017). This means that our model does not capture interannual variability due to freshwater inflow from rivers. We also omit tides in our model. River run-off and tides may contribute to interannual sea level variability but require a model setup with a higher horizontal resolution to be fully resolved.

Satellite Altimetry Data Used for Model Evaluation
We evaluated model results against sea level observations from the Level-4 gridded Ssalto/Duacs multimission satellite product ([1/4]°by [1/4]°) of AVISO (AVISO, 2019). This product is available for 1993 to a few months before present. It has been corrected for the IB effect and for the nodal cycle through standard tidal corrections. Furthermore, we assume that barystatic effects on interannual sea level variability are small compared to the effects of regional ocean dynamics, allowing for direct comparison with ROMS.
In addition, we evaluated our simulation of mean dynamic topography (MDT), which reflects the average strength of geostrophic currents. We used the CNES MDT CLS18 product (CNES, 2019), which is an estimate of the mean SSH above the geoid, and is given at a (1/8)°by (1/8)°resolution. It is based on altimetry data, in situ measurements of temperature and salinity, and the GOCO5S geoid model and covers the period 1993-2012.
Although the satellite altimetry data is given on [1/4]°and [1/8]°grids, respectively, its effective resolution is lower than our model resolution. Hence, small-scale features in the model results are difficult to evaluate. Additionally, land contamination leads to a lower accuracy of satellite altimetry in coastal regions (Andersen & Scharroo, 2011). Finally, satellite altimetry products are usually corrected for the IB effect, tides, and high-frequency winds (e.g., Ponte & Ray, 2002), leading to reduced accuracy in regions on the NWES where these corrections are difficult to model. Therefore, in section 3 we also compare our model results with GLORYS12v1.

ROMS Sensitivity Experiments
We decomposed sea level variability using six sensitivity experiments with different lateral and surface boundary conditions (Table 1). For the reference experiment, we varied all lateral and surface boundary conditions (Exp_AllVar). For the sensitivity experiments, we prescribed ROMS with an interannually varying annual cycle for the driver under investigation and with a repeated mean annual cycle for the others. We decomposed total interannual sea level variability (Exp_AllVar) into variability due to the lateral ocean boundary conditions, atmospheric forcing, and intrinsic variability (Exp_OceanVar, Exp_AtmosVar, and Exp_NoneVar). The influence of atmospheric forcing was further decomposed into the influences of SLP, winds, and buoyancy fluxes (Exp_SLPVar, Exp_WindVar, and Exp_BuoyVar). The SSH from GLORYS12v1 does not include the IB effect. Hence, when we include SLP in ROMS, the lateral boundary conditions from GLORYS12v1 are compensated using the PRESS_COMPENSATE option in ROMS.
All runs were initialized with the final year of Exp_NoneVar, that is, the sensitivity experiment in which both the atmospheric forcing and the lateral ocean boundary conditions are a repeated mean annual cycle. We calculated the mean annual cycle of the reanalysis data from the years 1993-1998, instead of the full period 1993-2018. This was done to minimize spin-up effects toward a changed climate state, as the mean state of Note. We used a climatological annual cycle for river run-off in every experiment.

Journal of Geophysical Research: Oceans
1993-1998 is more compatible with the boundary conditions in 1993. Additionally, we discarded the first two model years (1993 and 1994) to avoid any remaining spin-up effects. Therefore, in this manuscript we show model results for 1995-2018 (24 years). Spun-up sensitivity experiments show that results depend marginally on the period used to derive the climatological forcing.
Some of the bulk parameterizations of atmospheric forcing (section 2.2) depend on variables that vary in all sensitivity experiments. Thus, buoyancy fluxes may still vary interannually in sensitivity experiments other than Exp_BuoyVar. For example, the variability of SLP (Exp_SLPVar) affects air density and specific humidity, which in turn affect buoyancy fluxes. Since these feedbacks are ultimately caused by the interannual variability of the driver under investigation, we do not consider them part of Exp_BuoyVar.
Additionally, the response of the ocean to the various components of forcing at the lateral and surface boundaries depends on its state (temperature, salinity, and circulation). For example, the heat taken up by the ocean varies interannually in Exp_BuoyVar but is redistributed by the same wind-driven circulation each year. Hence, due to nonlinear feedbacks, the interannual sea level variability in Exp_BuoyVar may differ from the actual buoyancy-driven variability in Exp_AllVar. This applies to the other sensitivity experiments as well. However, we show in section 4 that such nonlinear feedbacks result in only small interannual sea level variability on the shelf.

Evaluation of Model Results
In this section, we evaluate simulated MDT against CNES MDT CLS18 (both as anomalies with respect to their area-weighted regional mean) and interannual sea level variability against the AVISO product. Additionally, we compare our simulations with the GLORYS12v1 reanalysis. We evaluate a variant of Exp_AllVar that excludes the effect of SLP. We compared the different fields quantitatively using pattern correlation coefficients (PCCs-the spatial equivalent of Pearson's correlation coefficient) and rootmean-square errors (RMSEs). We used a hypothesis test to determine statistical significance (p < 0.05, with p calculated using a two-sided t distribution).
ROMS MDT matches well with GLORYS12v1 MDT and CNES MDT CLS18 (Figures 2a, 2c, and 2e). It has PCCs of 0.95 and 0.93 and RMSEs of 3.76 and 4.34 cm with the other products, respectively. All three fields have a north to south gradient of low to high MDT anomalies. The MDT is higher than the regional mean on most of the shelf and largest in the southeastern North Sea. All fields show lower MDT in the Norwegian Trench and the northern North Sea. In ROMS and GLORYS12v1, the 0-cm MDT contour realistically follows the northern coast of Scotland and approximately traces the 80-m isobath into the North Sea.
Since GLORYS12v1 assimilates along-track satellite altimetry, its MDT is similar to the satellite altimetry data. However, in some areas, ROMS and GLORYS12v1 MDT match well but differ from CNES MDT CLS18. For example, ROMS and GLORYS12v1 differ only 1 cm near the coast of Scotland and Denmark, whereas CNES MDT CLS18 is up to 6 cm lower. This could be due to limitations of the satellite altimetry data (see section 2.3).
ROMS MDT differs from the other products near the Rockall Bank and in the southwest of the domain, so here geostrophic currents may not be well reproduced. Our model also has stronger MDT gradients across the shelf break, particularly at the Celtic and Armorican slopes. Such gradients indicate along-slope geostrophic currents, which have indeed been identified in this region with moorings and numerical models (Friocourt et al., 2008;Pingree & Le Cann, 1989) and may be affected by bathymetry smoothing. CNES MDT CLS18 shows no gradient in this area, but its resolution is likely insufficient to assess such a feature. Despite its data assimilation, GLORYS12v1 has a somewhat lower MDT along this part of the shelf break.
Other high-resolution regional ocean models with different boundary conditions also show, to various extents, reduced MDT along the Celtic and Armorican slopes (e.g., Hermans et al., 2020;Tinker et al., 2020).
Next, we compare the magnitude of interannual sea level variability in our model with the observations. To this end, we calculated the standard deviation (SD) of linearly detrended, annual mean sea level. Exp_AllVar reproduces the observed SDs well (Figures 2b, 2d, and 2f). It has PCCs of 0.80 and 0.78 and RMSEs of 0.49 and 0.53 cm with respect to GLORYS12v1 and AVISO. All three fields have large interannual sea level variability in deep waters, such as the Iceland Basin, the Rockall Trough, and the Bay of Biscay. However, the variability in ROMS (Figure 2b) is larger than in GLORYS12v1 and AVISO (Figures 2d and 2f) northwest of Portugal and generally smaller in the rest of the open ocean. Sea level variability is small near the shelf break and the Norwegian Trench, with all three fields showing a similar magnitude of 0.7-to 1.0-cm SD. This is encouraging, since the shelf break might play an important role for (decadal) sea level variability on the shelf .
The spatial patterns of interannual sea level variability also agree well on the shelf. As in GLORYS12v1 and AVISO, the largest variability in ROMS is found in the German Bight. The variability in the German Bight is somewhat higher in ROMS than in GLORYS12v1 and AVISO (up to 4.9-, 4.3-, and 4-cm SDs, respectively). In the central North Sea, ROMS is closer to AVISO, whereas AVISO has somewhat larger sea level variability in  Figure 5b). Intrinsic variability in ROMS does not have exactly the same spatial pattern as in the observations and is not expected to have the same phase. In the English Channel, the Irish Sea, and along parts of the shelf break, the correlations between ROMS and AVISO are lower than on the rest of the shelf, which could be due to limitations of the satellite altimetry data (see section 2.3). This is further supported by the pointwise correlation between GLORYS12v1 and AVISO, which is also reduced around the United Kingdom and parts of the shelf break (Figure 3b).
To further demonstrate the high correlations between sea level in ROMS and AVISO, we show the time series for a few sample TG locations at the coast (Figure 3c)  Finally, we compare the correlations of detrended sea level at Brest with detrended sea level at all other grid cells in ROMS and AVISO (Figures 4a and 4b), since the spatial correlation pattern in the eastern North Atlantic and on the NWES has been discussed in a number of previous studies (Calafat et al., 2012Chafik et al., 2019;Dangendorf et al., 2014;Frederikse, Riva, Slobbe, et al., 2016;Hogarth et al., 2020;Sturges & Douglas, 2011). The resulting patterns of ROMS and AVISO agree well within the model domain.
High correlations (>0.6) with sea level at Brest extend from northwestern Africa to Norway, with the highest correlations found primarily along the continental slope. In both ROMS and AVISO, sea level at Brest becomes progressively less correlated with sea level toward the southern and eastern North Sea. It is negatively correlated with sea level in the open ocean in the west and northwest of the domain, although this anticorrelation is mostly statistically insignificant. To highlight the unique character of the spatial coherence of coastal sea level variability, we compare the cross correlations of observed sea level along the continental slope and along a meridional transect in the open ocean (25°W) (Figures 4c and 4d). This shows that the correlation scales along the slope are several times longer than those in the open ocean, suggesting shelf-sea dynamics play an important role in driving spatial coherence on the NWES.
We conclude that, overall, interannual sea level variability in ROMS, GLORYS12v1, and AVISO is in good agreement (Figures 2-4). This indicates we can use our ROMS setup to investigate the different drivers of interannual sea level variability in the NWES region (section 4).

The Contributions of Different Drivers to Interannual Sea Level Variability
Using the sensitivity experiments in Table 1, we first isolate the contributions of intrinsic sea level variability, atmospheric forcing, and the lateral ocean boundaries (section 4.1). Then, we split up the contribution of atmospheric forcing into the effects of SLP, winds, and buoyancy fluxes (section 4.2). Finally, we show the relative contribution of different drivers to sea level at a few coastal locations (section 4.3) and relate these to the correlation between sea level at different locations (section 4.4).

Intrinsic Variability, Atmospheric Forcing, and Lateral Ocean Boundaries
In contrast to section 3, here we included the effect of SLP in Exp_AllVar (Figure 5a). The spatial patterns of sea level variability in Figures 2b and 5a are similar, but not exactly the same due to covariance. For example, including SLP leads to larger sea level variability in most locations but to slightly lower sea level variability in the central North Sea.
The interannual sea level variability in Exp_NoneVar (Figure 5b) has SDs of up to 2.7 cm in the deep ocean (Rockall Trough) but is negligible (<0.1-cm SD) on the shelf. Since all forcing in Exp_NoneVar is a repeated mean annual cycle, the resulting variability must be due to spontaneously generated nonlinear processes, such as mesoscale eddies. Therefore, we refer to this variability as intrinsic. A model with tides and The open ocean variability in Exp_NoneVar is lower than the intrinsic variability found in previous studies (e.g., Sérazin et al., 2015) and is negligible at the lateral boundaries. This is the result of the lateral boundary conditions, which were prescribed as a repeated annual cycle in Exp_NoneVar. The other sensitivity experiments also have intrinsic variability, so the sea level variability due to forcing at the lateral and surface boundaries is easier to identify on the shelf than in the open ocean.
Exp_AtmosVar (Figure 5c), in which only the atmospheric forcing varies interannually, shows that the atmospheric forcing is the dominant driver of sea level variability on the shelf. The resulting spatial pattern on the shelf compares well to Exp_AllVar (Figures 2b and 5a) and the observations (Figure 2f). Sea level variability in Exp_AtmosVar is highest in the German Bight (~5.6-cm SD) and is also relatively high around Scotland and Norway (~2.4-cm SD). Intrinsic sea level variability is present in the open ocean and is enhanced by the local effects of atmospheric variability. At the lateral boundaries, the interannual sea level variability is nudged toward the variability due to SLP, since a repeated annual cycle from GLORYS12v1 was used combined with the PRESS_COMPENSATE option (explained in section 2.4).
Exp_OceanVar (Figure 5d), in which only the lateral boundary conditions vary interannually, has only minor interannual sea level variability of around 0.5-cm SD on the shelf. The variability is mostly uniform on the shelf, as demonstrated by the high correlations between detrended sea level at Brest and elsewhere shown in supporting information Figure S1. By interannually varying the boundary conditions at each lateral boundary separately, we find that the sea level variability on the shelf is caused mostly by the variation of the western and southern boundaries of our model domain (see supporting information Figure S2). Variations in subtropical gyre strength may affect transport through the western and southern boundaries, but we did not find a significant correlation between sea level in Exp_OceanVar and the North Atlantic barotropic stream function in GLORYS12v1. In the open ocean, interannual sea level variability due to the varying lateral boundary conditions is much larger than on the shelf. Geostrophic currents along the continental slope can hinder the propagation of open ocean signals to the coast (Huthnance, 1981). In addition to the intrinsic sea level variability in the open ocean, the lateral ocean boundaries drive strong sea level variability along the western boundary, in the Iceland Basin (~4-cm SD) and east of the Faroe Islands (~2.6-cm SD).
Next, we calculated SDs for the linear sum of the detrended sea level time series in Exp_NoneVar, Exp_AtmosVar, and Exp_OceanVar (Figure 5e). We compared this to the variability in Exp_AllVar (Figure 5a), in which all forcing varies interannually simultaneously. The SDs for Exp_AllVar and the linear sum are similar on the shelf (Figures 5a and 5e), with differences well below 0.2 cm (Figure 5f). The linear sum also has the same phase, as indicated by the high pointwise correlation with Exp_AllVar (see supporting information Figure S3). This indicates that to first order, sea level variability due to the atmospheric and oceanic drivers combines linearly on the shelf. This also holds for the forcing at the separate lateral ocean boundaries (supporting information Figures S2 and S3). Since all sensitivity experiments include intrinsic sea level variability, differences between Exp_AllVar and the linear sum are larger in the open ocean (Figure 5f).
Exp_SLPVar (Figure 6b) shows the effect of SLP on interannual sea level variability. Except for the intrinsic variability in the open ocean, the resulting variability agrees with the variability of SLP in ERA5 (~1 cm sea level for 100-Pa SLP, see supporting information Figure S4). Thus, interannual sea level variability in Exp_SLPVar is mainly due to the IB effect. The resulting sea level variability is smallest in the south of the domain and in the southeastern North Sea (~1-cm SD) and increases toward the Icelandic Low northwest of the domain (~2.1-cm SD).
Exp_WindVar (Figure 6c)  Exp_BuoyVar (Figure 6d) results in spatially relatively smooth sea level variability that varies from 0.5-to 1.3-cm SD on the shelf. Apart from the intrinsic sea level variability, the magnitude of the variability in the open ocean is comparable to that on the shelf. Its magnitude tends toward 0 at the lateral boundaries because of the repeated annual cycle of GLORYS12v1. The interaction of buoyancy fluxes with river run-off and winds may explain the pattern of buoyancy-driven sea level variability in the German Bight.
Similar to section 4.1, we linearly summed the detrended sea level time series in Exp_SLPVar, Exp_WindVar, and Exp_BuoyVar and computed the resulting SDs (Figure 6e). The residual SDs relative to Exp_AtmosVar are less than 0.24 cm on the shelf (Figure 6f). Additionally, the pointwise correlation with the time series of Exp_AtmosVar is high (supporting information Figure S3). The residual SD is around 0.23 cm in the German Bight, which indicates that nonlinear feedbacks between buoyancy fluxes and winds, or buoyancy fluxes and river run-off, can have a small effect on interannual sea level variability in this area. The small effect of nonlinear feedbacks suggests that, to first order, the influences of different drivers of sea level variability on the shelf can be studied separately. This justifies the use of multiple (linear) regression methods in observation-based studies (e.g., Dangendorf et al., 2013Dangendorf et al., , 2014

Relative Magnitudes of Sea Level Variability at Coastal Locations
To investigate sea level variability in each sensitivity experiment, we display the time series for the sample coastal locations of Figure 3 (Figure 7, left panels). Additionally, we calculated the relative magnitudes of interannual sea level variability by dividing the SDs of detrended sea level in each sensitivity experiment by those in Exp_AllVar (Figure 7, right panels). These ratios do not necessarily add up to one for a given location, because the sea level variability in one sensitivity experiment may reinforce or interfere with that in another sensitivity experiment.
In most locations, SLP (red) drives the largest interannual sea level variability. This is especially the case for Brest, Santander, and Cascais (Figures 7f-7h, red). In Stavanger, Malin Head, and North Shields, both SLP (red) and wind (green) play a dominant role (Figures 7a-7c). Wind is the primary driver of interannual sea level variability at Cuxhaven and Den Helder (Figures 7d and 7e); there is large overlap between the time series of Exp_WindVar (green) and Exp_AllVar (black, dashed). The effects of winds at the western and eastern North Sea coasts clearly differ (Figures 7c-7e). Variability due to SLP, buoyancy fluxes, and the lateral boundary conditions is smaller than wind-driven variability at Cuxhaven and Den Helder but not

10.1029/2020JC016325
Journal of Geophysical Research: Oceans negligible. Using linear regression, Dangendorf et al. (2013) obtained similar results for the Cuxhaven TG record.
Since the interannual sea level variability in Exp_OceanVar is approximately uniform on the shelf (section 4.1), the relative importance of the forcing at the lateral boundaries depends on the pattern of sea level variability in Exp_AllVar. Hence, although sea level variability in Exp_OceanVar is relatively small (~0.5-cm SD, Figure 5d), the forcing at the lateral ocean boundaries is still a relatively important driver at Stavanger, North Shields, Santander, and Cascais. It results in SD ratios of 0.41, 0.33, 0.36, and 0.42, respectively (Figures 7a, 7c, 7g, and 7h, blue).
Especially Cascais is situated close to the lateral boundaries of our ROMS domain. Because of the nudging layer, sea level variability in Exp_WindVar and Exp_BuoyVar tends toward the repeated mean annual cycle of GLORYS12v1 at the lateral boundaries (Figures 6c and 6d). Therefore, the SD ratios of sensitivity experiments should be interpreted with caution at locations close to the lateral boundaries. Nevertheless, we expect that the drivers of sea level variability at Cascais, Santander, and Brest are similar, since sea level at these locations is highly correlated in both ROMS and the observations.

Correlation Between Sea Level at Coastal Locations
To relate the spatially varying effect of different drivers to the correlation between sea level at different locations, we calculated the correlation between detrended annual mean sea level in Exp_AllVar at several coastal locations in the region (Figure 8a, lower triangle). In section 4.2, we showed that wind variability is an important driver of interannual sea level variability in the North Sea ( Figure 6c). Consequently, sea level along the North Sea coast is highly correlated (e.g., at Wick, Aberdeen, North Shields, Esbjerg, Cuxhaven, Den Helder, and Oostende). Sea level at locations on the southwest of the NWES is also correlated (e.g., Brest, Newlyn, Santander, and Cascais), which is where SLP dominates (Figure 6b). Sea level at Malin Head is strongly influenced by both SLP and winds ( Figure 7b) and has high correlations with sea level at Bergen, Stavanger, and the northeast U.K. coastline. Sea level at Malin Head has reduced correlations with sea level at locations in the North Sea, since these are predominantly controlled by wind variability. Correlations between sea level at locations in the North Sea and the English Channel, and in the North Sea and on the southwestern NWES, are low or even negative (e.g., Santander and Brest).
The correlation pattern for Exp_AllVar generally agrees with that for PSMSL (Holgate et al., 2013;PSMSL,2020) TGs (Figure 8a, upper triangle), but TG correlations are generally lower, possibly due to data gaps and the influence of local processes. We excluded TGs with less than 18 years of data (75%, gray) and corrected for the nodal cycle following Woodworth (2012), neglecting self-attraction and loading effects. Both ROMS and the TGs show that sea level on the southwestern part of the NWES has low or negative correlations with sea level in the southern and eastern North Sea. This can be explained by the sensitivity of shelf sea level to the NAO (Chen et al., 2014;Tsimplis & Shaw, 2008;Tsimplis et al., 2005;Wakelin et al., 2003;Woolf et al., 2003). For example, a high SLP anomaly at Santander depresses sea level locally. This coincides with stronger southwesterly winds over the North Sea, which induce higher sea level at the North Sea coast through Ekman transport (Tsimplis & Shaw, 2008). This mechanism becomes apparent from Exp_SLPVar and Exp_WindVar (Figures 6b and 6c).
Complementing Figure 8a, we repeat Figure 4a, now incorporating the effect of SLP (Figure 8b). This leads to increased and statistically significant correlations between sea level at Brest and in the open ocean, because SLP has long correlation scales and the coastal region and open ocean respond to SLP similarly. As discussed earlier, SLP also enhances the contrast between sea level variability on the southwestern part of the NWES and in the North Sea. As a consequence, the coherence of sea level variability along the shelf break extends less far north in Figure 8b than in Figure 4a.
To investigate the role of LSWs inside our domain in driving the coherence of sea level on the NWES, we integrated ERA5 winds along the 400-m isobath from the southern lateral boundary (36°N) to the respective latitudes of Stavanger, Malin Head, Brest, Santander, and Cascais. The LSW integrals are significantly correlated with detrended annual mean sea level in Exp_AllVar (excluding the IB effect) at these locations, which primarily reflect high correlations with sea level in Exp_WindVar (supporting information Table S1). The correlations of sea level with LSWs are higher than those with local winds, supporting the premise of signals propagating along the boundary (Calafat et al., 2012;Sturges & Douglas, 2011). The

Journal of Geophysical Research: Oceans
LSW integral up to Stavanger has relatively low correlations with sea level at North Shields, Cuxhaven, and Den Helder. In line with Figure 8b, including SLP affects the correlation between integrated LSWs and annual mean sea level (supporting information Table S1), reflecting the important role of the IB effect at most locations.

Discussion
Our sensitivity experiments showed that atmospheric forcing is the main driver of interannual sea level variability on the NWES. This agrees with previous studies for different time periods (e.g., Dangendorf et al., 2014;Tinker et al., 2020;Wakelin et al., 2003). Additionally, we found that ocean-driven variability is small and coherent on the shelf and combines with atmospheric-driven variability nearly linearly. Similar spatial patterns of sea level variability were obtained by Tinker et al. (2020) with a different regional ocean model, indicating some robustness to model design. Because Tinker et al. (2020) derived boundary conditions from a 200-year control run of HadGEM3 GC3.0 (Williams et al., 2018), the magnitude of sea level variability in their atmospheric and ocean-driven simulations differs from that in Exp_AtmosVar and Exp_OceanVar (Figures 5c and 5d). Since we used reanalysis data as boundary conditions instead, our results can be compared directly with satellite altimetry but may not be fully representative of (multi)decadal sea level variability.
We have separated the atmospheric forcing into SLP, winds, and buoyancy fluxes for the entire NWES. This complements the analysis of Dangendorf et al. (2014), which focused on TGs in the North Sea. Dangendorf et al. (2014) found that linear regression with winds explains around 50-85% of observed variability from Oostende (Belgium) to Hirtshals (Denmark) and with SLP around 50-65% from Tregde to Bergen (Norway) and Aberdeen to Lerwick (UK). Exp_WindVar confirms that wind-driven interannual sea level variability is largest in the southern and eastern North Sea. Our results also show that wind-driven variability is large around Norway, Scotland, and north of Ireland and small on the southwest of the NWES and along the shelf break. In line with Dangendorf et al. (2014), we showed that SLP is the dominant driver of interannual sea level variability along the northeast U.K. coastline. Our sensitivity experiments show that this is also the case on the southwestern part of the NWES. Much of the interannual variability at North Sea TGs is explained with barotropic models Piecuch et al., 2019). The baroclinic response to winds could be further investigated by comparing our results to the results of a barotropic setup under the same atmospheric forcing.
Since the NWES seas are shallow, local steric sea level variability is limited (Roberts et al., 2016;Tinker et al., 2020). We showed that buoyancy fluxes result in interannual sea level variability with SDs varying from 0.5-1.3 cm on the shelf (Exp_BuoyVar; Figure 6d). At the sample locations in Figure 7, the relative magnitude of buoyancy-driven variability roughly agrees with the finding of Tinker et al. (2020) that the local steric component explains up to 35% of the total variability on the shelf. A large fraction of the steric sea level variability on the NWES in the global model of Roberts et al. (2016) was indeed explained by buoyancy fluxes.
The coherence of sea level along the eastern boundary of the North Atlantic and in the North Sea was previously attributed to wind-driven CTWs (Calafat et al., 2012Dangendorf et al., 2014;Sturges & Douglas, 2011). The long correlation scales along the shelf break compared to in the deep ocean (Figures 4c and 4d) suggest that boundary dynamics are indeed at play. This is further supported by the significant correlations of integrated LSWs with simulated coastal sea level (supporting information Table S1), which are similar to the correlations with TG records (Calafat et al., 2012Dangendorf et al., 2014). The correlations between integrated LSWs and sea level in Exp_WindVar confirm that this results from wind forcing, and the reduced correlations between local winds and sea level suggest that LSWs drive an important part of the variability in Exp_WindVar. The correlation between integrated LSWs and sea level in the North Sea is low, mirroring the correlation pattern of a common mode of British TGs with satellite altimetry in (Hogarth et al., 2020) and indicating the dominance of local winds over the North Sea (Figure 6c).
Since wind-driven interannual sea level variability is relatively small at the shelf break and on the southwestern NWES (Exp_WindVar, Figure 6c), the coherence of sea level there could partially be explained by the relatively uniform sea level variability due to the forcing at the lateral ocean boundaries and due to buoyancy fluxes (Exp_OceanVar and Exp_BuoyVar, Figures 5d and 6d). The correlation scales in Figure 4 and the correlation between sea level and integrated LSWs strongly suggest that wind-driven CTWs influence sea level on the NWES but do not necessarily rule out an effect of open ocean steric anomalies (Chafik et al., 2019), which may be represented in Exp_OceanVar. The variability driven by the southern boundary conditions may also reflect northward propagating CTWs triggered south of the domain (Calafat et al., 2012;Fukumori et al., 2015). However, quantifying this remote contribution accurately would require enlarging the model domain and increasing the simulation period.

Conclusions
We have used a high-resolution ([1/8]°by [1/8]°) configuration of ROMS for the NWES to identify the contribution of different drivers of interannual sea level variability during 1995-2018. The interannual sea level variability in our model matches well with that observed with satellite altimetry. Based on sensitivity experiments in which we isolated the effect of atmospheric and oceanic drivers, we conclude that atmospheric variability is the main driver of interannual sea level variability, with the largest SDs in the German Bight (up to 5.6 cm). On the other hand, the forcing at the lateral ocean boundaries causes coherent sea level variability on the shelf, with an SD of only 0.5 cm. The ocean-driven variability is mainly caused by the variation of the western and southern boundary conditions, but its driving physical mechanism requires further investigation. In our model, intrinsic sea level variability has an SD of up to 2.7 cm in the open ocean but is negligible on the shelf.
Our results show that, to first order, the interannual sea level variability induced by variability of SLP, winds, buoyancy fluxes, and the lateral ocean boundaries combines linearly on the NWES. Thus, the contributions of these drivers can be studied separately, which supports the use of multiple linear regression analysis. Wind variability has a large influence on sea level variability in the North Sea and around Scotland, Norway, and north of Ireland. It is the dominant driver in the southern and eastern North Sea, whereas the IB effect dominates at other coastal locations shown in this study. Wind also drives large mesoscale variability in the open ocean, but wind-driven sea level variability is suppressed along the shelf break. At the sample locations in Figure 7, buoyancy-driven sea level variability can have a relative magnitude of up to 55% (North Shields). Sea level variability driven by the lateral ocean boundary conditions can have a relative magnitude of up to 42% (Cascais).
Sea level variability on the southwestern part of the NWES is predominantly driven by the IB effect. As a consequence, detrended annual mean sea level on the southwestern part of the shelf has low or even negative correlations with that in the southern and eastern North Sea, where sea level variability is predominantly driven by local winds. This agrees with the sensitivity of shelf sea level to the NAO (e.g., Tsimplis & Shaw, 2008;Wakelin et al., 2003). The significantly positive correlations between wind-driven sea level variability on the NWES and winds integrated along the continental slope in our model indicate that LSWs and the consequent propagation of signals along the eastern boundary of the North Atlantic are partly responsible for the coherence of sea level variability on the NWES.
To conclude, our results improve the understanding of the drivers of interannual sea level variability on the NWES, which vary in importance depending on location. Explicitly decomposing the contribution of different drivers using ROMS increases confidence in previous observation-based studies and is useful for sea level budget studies. The results of our sensitivity experiments can be used to remove known variability from observations to estimate sea level rise or time of emergence with higher accuracy. This will help coastal decision makers to monitor and detect early warning signals more reliably.