Dynamical downscaling of unforced interannual sea-level variability in the North-West European shelf seas

Abstract

Variability of Sea-Surface Height (SSH) from ocean dynamic processes is an important component of sea-level change. In this study we dynamically downscale a present-day control simulation of a climate model to replicate sea-level variability in the Northwest European shelf seas. The simulation can reproduce many characteristics of sea-level variability exhibited in tide gauge and satellite altimeter observations. We examine the roles of lateral ocean boundary conditions and surface atmospheric forcings in determining the sea-level variability in the model interior using sensitivity experiments. Variability in the oceanic boundary conditions leads to uniform sea-level variations across the shelf. Atmospheric variability leads to spatial SSH variability with a greater mean amplitude. We separate the SSH variability into a uniform loading term (change in shelf volume with no change in distribution), and a spatial redistribution term (with no volume change). The shelf loading variance accounted for 80% of the shelf mean total variance, but this drops to ~ 60% around Scotland and in the southeast North Sea. We analyse our modelled variability to provide a useful context to coastal planners and managers. Our 200-year simulation allows the distribution of the unforced trends (over 4–21 year) of sea-level changes to be quantified. We found that the 95th percentile change over a 4-year period can lead to coastal sea-level changes of ~ 58 mm, which must be considered when using smooth sea level projections. We also found that simulated coastal SSH variations have long correlation length-scales, suggesting that observations of interannual sea-level variability from tide gauges are typically representative of > 200 km of the adjacent coast. This helps guide the use of tide gauge variability estimates.

1 Introduction

Relative sea-level change is one of the most important aspects of a changing climate. In addition to an anthropogenic driven trend in present and future sea level (Church et al. 2013; Marcos and Amores 2014; Slangen et al. 2014, 2016) variability may show itself as short term accelerations in sea-level rise (Calafat and Chambers 2013; Haigh et al. 2014; Marcos et al. 2017). Recent sea-level projections developed for the UK (Palmer et al. 2018) suggested that for a given tide gauge, variability will dominate over the emission scenario uncertainty and model structural uncertainty for the next decade (e.g., Newlyn, UK, Fig. 1) and will remain an important component throughout the twenty-first century.

Fig. 1

Adapted from Palmer et al. (2018)

Regional sea-level projections for Newlyn (UK) to illustrate the sea-level variability relative to the projected change and uncertainty. All time-series are plotted relative to a baseline period of 1981–2000. Coloured lines indicate the central estimates according to the figure legend. Shaded regions represent the projection range for the corresponding relative concentration pathway (RCP, left panel). The fraction of sea-level rise uncertainty from: sea-level variability (yellow); climate change scenario (green); and model uncertainty (blue), following Hawkins and Sutton (2011) (right panel).

Sea-surface height (SSH) variability has been studied with tide gauge and other in situ measurements (e.g., Tsimplis and Shaw 2008; Wahl et al. 2013; Frederikse et al. 2016), satellite observations (e.g., Hakkinen and Rhines 2004; Chafik et al. 2019) and models (e.g., Chen et al. 2014; Roberts et al. 2016). There are several components of SSH that can lead to variability, including steric sea-level change, mass convergence (manometric sea-level change, Gregory et al. 2019) and the inverted barometer effect (Stammer and Hüttemann 2008). Each of these terms can vary on different time and space scales and have a different response to modes of climate variability—see Roberts et al. (2016) for a review. Furthermore, shelf exchange processes can mediate how tightly coupled these components are between the ocean and the adjacent shelf seas (e.g., Landerer et al. 2007; Bingham and Hughes 2012; Chafik et al. 2019).

An understanding of the processes that lead to variability may allow prediction of SSH variability (e.g., Roberts et al. 2016; Sonnewald et al. 2018). Advanced warning of extreme sea level events is an invaluable tool for coastal communities, allowing the implementation of management policies and strategies to minimise loss of life and infrastructure damage (Miles et al. 2014). There have been several studies that have investigated the predictability of seasonal-to-interannual sea-level variability, based on statistical (e.g., Xue and Leetmaa 2000; Chowdhury et al. 2007, 2014; Sonnewald et al. 2018) and dynamic models (e.g., Miles et al. 2014; McIntosh et al. 2015; Roberts et al. 2016; Widlansky et al. 2017; Sonnewald et al. 2018). There tends to be more predictability in the tropics (where the SSH varibility is largely steric in nature) than in the extratropics and shallow marginal seas (such as the northwest European shelf seas) where wind-driven variability is more important (Roberts et al. 2016). The atmosphere has less memory than the ocean (Roberts et al. 2016), and so SSH variability associated with wind-driven, rapid barotropic adjustment is less predictable (Miles et al. 2014; Roberts et al. 2016). For example, Häkkinen (2004) showed that in the eastern Atlantic (around Europe) local and remote wind stress forcing plays an important role in SSH variability. Miles et al. (2014) showed that correlations between a corrected global reanalysis and altimeters are ~ 1 in the tropics, but < 0.5 around Europe. Seasonal predictability of the North West European Shelf Seas (NWS) marine environment (including SSH) is in its infancy (Tinker et al. 2018) but is an active field of research.

Europe, and particularly the North Sea coast, has several long tide gauge records that allow an observation-based assessment of sea-level variability. Dangendorf et al. (2014) examined interannual to decadal variability in 22 European long tide gauges records (some starting in the late nineteenth century) and considered different frequencies and their possible drivers. They found decadal-scale variability similar among all tide gauges, while the higher frequency variability (with timescales up to a couple of years) varied between tide gauges and was related to atmospheric drivers. The high frequency variability of tide gauges in the North Sea was predominately atmosphere driven, with a region dominated by the inverse barometer effect spanning from the English Channel diagonally towards Norway, and a wind stress driven region in the South-Eastern North Sea. The lower frequency variability was partly explained by the remote steric effect originating from adjacent to the shallow North Sea (in the Norwegian Trench). On the decadal timescale, wind stress-driven coastally trapped waves (which can travel large distances along the shelf break) have been hypothesised to make an important remote contribution (Sturges and Douglas 2011; Calafat et al. 2012). By characterising and removing these terms it is possible to calculate the long-term trend and acceleration in the tide gauge records more accurately. Frederikse et al. (2016) were able to close the sea-level budget for (a virtual mean of) the tide gauges around the North Sea, for the trend and the variability. On the interannual to decadal timescale, a number of studies have shown the importance of climate drivers such as the North Atlantic Oscillation (NAO) (e.g., Wakelin et al. 2003; Tsimplis et al. 2006; Chen et al. 2014) and state of the sub polar gyre (Chafik et al. 2019) in driving NWS sea-level variance. For example, Wakelin et al. (2003) compared the winter mean sea level and NAO and showed the southeast North Sea had correlations > 0.7 with a sensitivity of up to 96 mm/unit NAO index.

The North West European Shelf Seas (NWS) are the broad shallow seas to the north-west of Europe. They are bounded by several populous countries, including the UK, Ireland, France, The Netherlands, Germany and Norway, and are of significant economic, cultural and environmental importance. The NWS are quasi-isolated from the adjacent North-Eastern Atlantic by the steep gradient of the shelf slope, and the shelf slope current (Wakelin et al. 2009), and so the oceanographic conditions on and off the shelf are very different and evolve separately. They have particularly energetic tides (Pugh 1987) that determine whether regions seasonally stratify or remain fully mixed throughout the year (Simpson and Bowers 1981; Elliott and Clarke 1991). Current (e.g., CMIP5 and CMIP6, the 5th and 6th phase of the Coupled Model Intercomparison Project, Taylor et al. 2012; Eyring et al. 2016) Global Climate Models (GCMs) that form the basis of most process-based global and regional sea-level projections (e.g., Slangen et al. 2012, 2014; Church et al. 2013; Kopp et al. 2014; Cannaby et al. 2016) do not include dynamic tides, and so the NWS are typically poorly represented in terms of temperature and salinity (e.g., Mathis et al. 2013), and sea level (e.g., Hermans et al. 2020b). For these reasons, we do not analyse GCMs SSH variability directly (although we do compare the performance of the NWS SSH interannual variability with and without downscaling, in the additional materials).

Driving shelf-seas models with output from GCMs (dynamic downscaling) allows a more realistic simulation of the NWS through improved horizontal and vertical resolution and the inclusion of important processes typically neglected in CMIP6 GCMs (such as tides). This approach is well established for projections of temperature and salinity (e.g., Ådlandsvik 2008; Holt et al. 2010; Tinker et al. 2016) but is less common for studies of sea level (e.g., Olbert et al. 2012; Chen et al. 2014; Mathis and Pohlmann 2014; Hermans et al. 2020b). Here we use a regional shelf seas model (NEMO Shelf Coastal Ocean model version 6, CO6, O’Dea et al. 2017) to investigate unforced year-to-year variability on the NWS. In a related study, Hermans et al. (2020b) used our model set up to downscale a pair of CMIP5 models to investigate projections of sea level. We downscale a “present day” (conditions representative of the year 2000) control simulation of the Met Office Hadley Centre CMIP6 global coupled model HadGEM3 GC3.0 (Williams et al. 2018). Climate control simulations provide insights into unforced climate variability by providing longer time-series than the available observational records. We evaluate the last 200 years of the simulation against observed sea-level variability from satellite and tide gauge records. We focus on interannual-to-decadal timescales, but note that there is important sea level variability on longer timescales (e.g., decadal-centennial and longer), associated with processes such as the Atlantic Multidecadal Variability (McCarthy et al. 2015) and barystatic sea-level change (Frederikse et al. 2016).

After showing how well our modelling system can reproduce interannual variability on the NWS, we aim to investigate the behaviour of the modelled SSH variability within the NWS and explore the relative roles of the variability associated with the ocean boundary conditions, and the atmosphere surface forcings. We do this by separating the SSH into a spatially coherent shelf loading term, and a redistribution term. We investigate the relative contributions of these terms in different locations, and on different time scales. As this work was funded by, and fed into, the marine component (Palmer et al. 2018) of the United Kingdom’s Climate Projections of 2018 (UKCP18), an important secondary aim of this study is to describe the modelled variability in a practical manner, to provide context to sea-level projections and observation-based estimates of sea level variability. We do this by assessing how much sea-level change could be expected due to unforced variability within a given period, and how representative observation-based estimates of sea level variability are of the adjacent coast.

2 Models and methods

This study makes use of two independent modelling systems in a “nested” configuration (Fig. 2, Table 1). The first is the GCM HadGEM3, version GC3.0 (GC3.0, Williams et al. 2018, which is essentially the same physical model (i.e., atmosphere, ocean and sea-ice) submitted to CMIP6 by the Met Office). This model (Table 1) represents the state-of-the-art in coupled climate modelling. Output from GC3.0 is used to provide surface fluxes and lateral boundary conditions for a regional model of the NWS (Fig. 2). This second model system is NEMO 3.6 Coastal Ocean model version 6 (CO6, O’Dea et al. 2017) and represents the state-of-the-art in coastal ocean modelling and is the basis of the Copernicus Marine Environmental Monitoring Service (CMEMS) NWS products. Details of these modelling systems and the specification of boundary conditions are provided in the following sections.

Fig. 2

GC3.0 and CO6 model grids, showing model coastline (black), depth contours (at 50, 100, 200, 500, 1000 and 2000 m, grey), and every 5th ocean grid box (blue dots). a the GC3.0 model grid, also showing every 5th atmosphere grid box (red dots). b CO6 model domain, also showing the extent of the NWS (red)—it is (approximately) delineated by the 200 m isobath, but excludes the Armorican shelf, the Skagerrak/Kattegat and the Norwegian Trench

Table 1 Model description

2.1 Model description

In this section we provide a brief overview of the HadGEM3 GC3.0 and NEMO 3.6 CO6 models in terms of their key features and refer the reader to the literature for further information.

2.1.1 The HadGEM3 GC3.0 coupled climate model

HadGEM3 GC3.0 is the Met Office Global Coupled model version 3.0 (Williams et al. 2018). The atmosphere component uses the Met Office Unified Model with GA7.0 settings (Walters et al. 2019) at N216 horizontal resolution (~ 60 km) and 85 vertical levels. The Boussinesq ocean component uses the NEMO model with GO6.0 settings (Storkey et al. 2018) at ORCA025 resolution (~ 1/4° on a tri-polar grid) and 75 vertical levels. The model has a non-linear free surface (see appendices) but does not simulate tides or the inverse barometer effect. The sea ice component uses the CICE model (Hunke et al. 2015) with GSI8.0 science (Ridley et al. 2018) and uses the same tri-polar grid as NEMO.

For this study we use GC3.0 to run a present-day control simulation. It is initialised from EN3 climatology (2004–2008) of temperatures and salinity (Ingleby and Huddleston 2007) and run for 270 years (model years 1980–2250). During this time greenhouse gas concentrations, ozone concentrations and aerosol emissions are kept constant (or with a repeating annual cycle) at year 2000 levels, and hence this simulation represents the near-present day for the duration of the model run (e.g., model years 1980 and 2250 both represent conditions consistent with the year 2000).

This present-day control simulation has been assessed (for the years 50–100) by Williams et al. (2018). The GC3 global mean thermosteric sea level has a linear trend of 0.599 mm/year, which is much lower than the present day or recent past value of global mean sea level change, but suggests the deep ocean is still spinning up. The GC3.0 Atlantic Meridional Overturning Circulation (AMOC) is 18.6 Sv (± 1 standard deviation = 1.0 Sv) at 26° N (1000 m) which is comparable to the observed value of 17.0 Sv ± 4.4 Sv (Frajka-Williams et al. 2019), and leads to a Northward heat transport within the observation uncertainty of Ganachaud and Wunsch (2003). We further assess the climatological precipitation and wind stress biases over the region of the CO6 model for each seasonal mean in the GC3.0 driving model and for comparison, the previous CMIP5 generation Met Office climate model. The root mean square (RMS) biases over the CO6 region are calculated by comparing the model climatologies with GPCP precipitation (Adler et al. 2003) and SCOW wind stresses (Risien and Chelton 2008) (see Table 8 in additional materials). GC3.0 reduces wind stress biases over this region making this a superior model for driving CO6. Precipitation biases are better in June to August but worse in December to February; and with biases less than 1 mm day⁻¹, GC3.0 is still sufficiently accurate to provide surface forcings for the CO6 model. The NAO (a primary mode of seasonal to decadal variability in the NWS region) has improved in the GC3.0 model (with respect to the previous model) as shown by Scaife et al. (2011) where increasing the ocean resolution to ORCA025 (as is used in GC3.0) improved temperatures in the North Atlantic which improved both the North Atlantic blocking and seasonal forecasts of the NAO (Scaife et al. 2014). Given this, and the agreement of the modelled and observed sea-level variability (below) we consider GC3.0 to be fit of purpose as a driving GCM.

2.1.2 The 7 km European shelf seas model

NEMO Coastal Ocean model version 6 (CO6) implementation (O’Dea et al. 2017) is a primitive equation, Boussinesq, 3D baroclinic model, with a non-linear free surface (see Sect. 11.1). It is run on a regional ~ 7 km grid extending from 40° 4′ N 19° W to 65° N 13° E, with 50 hybrid terrain following vertical levels (s-levels, Siddorn and Furner 2013). This resolution is insufficient to resolve the internal (baroclinic) Rossby Radius on the shelf (which is of the order 4 km) but resolves the external (barotropic) Rossby Radius (∼ 200 km). Ideally the model would be of sufficient resolution to resolve both the internal and external radii, i.e., have a resolution of the order < 2 km (O’Dea et al. 2012), however, the computation expense of such resolution is impractical for climate integrations. Both tides and the inverse barometer effect are modelled directly and are also applied to the lateral ocean boundary conditions.

CO6 is a well-established and evaluated model, with a wide range of uses. It is used as a research model, as the basis of the Met Office operational 6-day NWS forecasts (and delivered to CMEMS (https://marine.copernicus.eu/services-portfolio/access-to-products/), Tonani et al. 2019) and 26-year NWS reanalyses (also delivered to CMEMS, Renshaw et al. 2019).

2.1.3 The Boussinesq approximation

Both the NEMO ocean component of HadGEM3 GC3.0 and CO6 make the Boussinesq approximation. In the model equations the in situ density is replaced by a constant reference density in all but the vertical momentum equations and the equations of state (Gill 1983). As a result, these models conserve volume rather than mass, and so do not model all components of the steric expansion prognostically.

The steric component of sea level change is the combination of three effects: the local steric effect; the global-mean steric effect; and the non-Boussinesq steric effect (Griffies and Greatbatch 2012). The local steric effect contributes to sea level through the Eulerian time derivative of depth-integrated local in situ density (how density changes in a fixed location), and is included in both Boussinesq and non-Boussinesq models (Griffies and Greatbatch. 2012). The global steric effect gives rise to global mean sea-level change through changes in global mean density. The non-Boussinesq steric effect arises from the depth-integrated Lagrangian time derivative of depth-integrated in situ density (how the depth-mean density changes as you follow a water parcel). The sea level in non-Boussinesq models includes all three of these steric effects. In Boussinesq models, non-Boussinesq steric effects and the global steric effect are not modelled prognostically, but can be diagnosed a-posteriori (Greatbatch 1994).

In this study, our models therefore include the local steric effect, but do not include the global mean and non- Boussinesq steric effect. We do not correct for these effects, as we expect that they do not greatly impact regional sea-level variability on frequencies larger than monthly (Griffies and Greatbatch. 2012). Since the local steric effect arises from the time tendency of depth-integrated density anomalies, a low-density anomaly in the open ocean will lead to a greater local steric effect than in the adjacent shelf seas. This will lead to a horizontal pressure gradient and barotropic adjustment, driving water from the open ocean into the shelf seas (Bingham and Hughes 2012).

2.1.4 Shelf seas model forcings

We follow essentially the same dynamic downscaling approach as described by Tinker et al. (2015, 2016) and Hermans et al. (2020b), using model output from a GCM to generate model forcings for the shelf seas model in a one-way nested configuration. There are broadly four types of model forcings: atmospheric surface boundary forcings; oceanic lateral boundary forcings; river forcings; and Baltic Sea exchange forcings (see Table 9 in the additional materials for details of how these are implemented). We generate oceanic and atmospheric forcings based on time-evolving fields simulated by GC3.0 and use climatological forcings for river run-off and the Baltic Sea. An additional simulation is run using interannually varying river run-off modelled by GC3.0, but this has a little effect on sea-level variability. The use of climatological boundary conditions for the exchange between the Baltic and North Sea is well established for NWS downscaling studies (Ådlandsvik 2008; Holt et al. 2010; Chen et al. 2014; Mathis and Pohlmann 2014; Tinker et al. 2016). Furthermore, Hermans et al. (2020a submitted) show that variability in the Baltic boundary conditions has a small effect on the NWS SSH variability. Given the cyclonic circulation of the North Sea, the impact of (and reduced variability associated with) the climatology for the Baltic exchange is limited to the Norwegian Trench (e.g., Holt et al. 2010; Tinker et al. 2016; Hermans et al. 2020a submitted).

Surface fluxes of heat, freshwater and momentum are taken directly from the GC3.0 atmosphere and interpolated onto the CO6 grid using a Gaussian interpolation scheme. GC3.0 provides instantaneous hourly wind and pressure data, 3-h mean Sea Surface Temperature (SST from GC3.0’s ocean), and 6-h mean heat and freshwater fluxes. Daily-mean ocean temperature and salinity (3D fields), surface elevation anomaly (relative to the instantaneous global mean) and barotropic currents (2D fields) are taken from the GC3.0 ocean and interpolated onto the CO6 grid at the lateral boundaries of the model domain. When creating the CO6 surface elevation anomaly boundary conditions, we do not add the GC3.0 global mean thermosteric sea level (“zostoga” in CMIP parlance) to the GC3.0 dynamic sea level (zos). Exchange with the Baltic Sea is treated as an oceanic lateral boundary condition, with T, S and barotropic currents specified from a model-based climatology (O’Dea et al. 2017). Surface elevation is not prescribed in the Baltic Sea exchange boundary conditions to avoid inconsistencies with the SSH specified at the Atlantic open boundary. We use a similar climatological river dataset to Graham et al. (2018), providing daily mean river volume fluxes. Nineteen tidal constituents, calculated from a tidal model of the North-East Atlantic (Flather 1981), are specified for the boundary depth mean velocities and sea surface elevation. As the CO6 domain region covers a significant area, the equilibrium tide (tidal potential) is also specified (O’Dea et al. 2012). We use a Flather radiation boundary condition (Flather 1976) to allow information to freely propagate in and out of the domain, while maintaining the realistic simulation of NWS tides. Initial conditions (of 3D temperature and salinity) are interpolated from GC3.0 with a nearest-neighbour interpolation. Wakelin et al. (2009) have shown that initial conditions can influence the NWS model solution for up to 7 years. We avoid this with a spin-up period (see below).

2.1.5 Experimental design

The GC3.0 present-day control simulation begins in 1980 and runs forward with fixed year 2000 greenhouse gas concentrations. The model runs until 2250, and the full control simulation is downscaled with CO6. The downscaled simulation exhibits drift in SST, SSS (sea-surface salinity), and SSH within the first 70 years, associated with the GC3.0 spin-up (Table 11 in additional materials gives the drift in the first 70 years, and in the subsequent 200 years of the downscaled simulation). In this study we therefore consider the first 70 years as spin-up and focus our analysis on the last 200 years (model year 2050–2250). This simulation is referred to as “Ctrl”.

A pair of companion sensitivity tests are run to compare the importance of the atmospheric and the ocean variability for the NWS. We create a 30-year mean seasonal cycle as a climatology of the atmospheric and oceanic forcings (2041–2070). One simulation (starting from year 2050 from the main simulation, Ctrl) uses the climatological atmosphere and the original (i.e., time evolving) ocean forcings to isolate the impact of the ocean boundary variability on the NWS SSH variability (CtrlOcV). A second sensitivity simulation uses climatological ocean forcing, with the original atmosphere forcings to isolate the impact of atmospheric variability (CtrlAtV).

2.2 Data

We use observations and analysis products to evaluate the model run. As our simulations are driven by a climate model, the phase of climate variability is not expected to match that of the observations. Our model evaluation therefore focuses on a comparison of statistics of model simulated and observed SSH variance for the tide gauge and satellite altimetry products. We focus model evaluation and analysis on the last 200 years of the simulation. We also include a broader model evaluation, in terms of temperature and salinity, in the additional materials section.

2.2.1 Tide gauges from permanent service for mean sea level

We evaluate the model SSH variability with tide gauge data from the Permanent Service for Mean Sea Level (PSMSL, Holgate et al. (2013); Data retrieved 20th of September 2019). We consider the modelled SSH from the model grid box nearest to the tide gauge location to be the model-simulated tide gauge. We use PSMSL monthly revised local reference (RLR) data to account for changing baselines and we reject data with quality issues (calculations for mean tide level, suspect data, etc.). Annual mean time-series are created from monthly mean anomalies (where the climatological season cycle has been removed) where there is data from at least 11 months. We use tide gauges within the NWS, with at least 15 years of data that meet this requirement. As many tide gauge records have breaks, we use the longest continuous section for spectral analysis—this is typically much shorter than the full record (see Table 2). As we are interested in interannual to decadal variability, we remove the linear trend from the data. Furthermore, when spectral analysis is not performed (i.e., Figs. 3 and 4), we high pass filter the data with a 32-year threshold—this is considered longer than any period of interest. We use standard deviations calculated from tide gauges with at least 15 years of data (after quality control and removal of years with less than 11 months of data) for our spatial variability comparison (Fig. 4). For our spectral analysis, we restrict our analysis to six tide gauges with continuous records greater than 50 years.

Table 2 Summary statistics of the observed and model simulated 32-year high pass filtered SSH anomaly (mean has been subtracted) from Fig. 3, including the median sample standard deviation, skewness and kurtosis

Fig. 3

Comparison of selected (32-year high pass filtered) observed (green) and model simulated (black) SSH anomaly (time-series mean removed) at tide gauges around Europe (with the bracketed letter in the title corresponding to the location on Fig. 4). The observed tide gauge record is plotted against time. The model-simulated SSH anomaly is appended afterwards. Horizontal green and black lines show the ± 1.96 standard deviation above and below the mean (grey line, 0 mm) for the observed and modelled datasets respectively. On the left of each time-series is a representation of the two distributions (estimated with a kernel-density estimate with Gaussian kernels and a bandwidth of 0.3) with the median value plotted as horizontal bold line. The distributions and median line illustrate how symmetrical the distribution is. Median, sample standard deviation, skewness and kurtosis are given for each tide gauge in Table 2

Fig. 4

Spatial comparison of variability of the quality controlled, observed and model-simulated tide gauges. Both time-series are linearly detrended and 32-year high pass filtered. Only locations with 15 years of data are included. The model data are not sub-sampled to match the length of the tide gauge record. Left: tide gauge variability (annual mean sample standard deviation) are illustrated as bar graphs at the tide gauge location. The left bar (green) is the observed value while the model-simulated value is the right (black) bar. Right: scatter plot showing the linear relationship (blue line) between these observed and model-simulated tide gauge data, with a 1:1 line (grey) and statistics of this relationship given [Correlation (r), relative standard deviation (rsd), and the coefficients of the linear fit (y = mx + c)]. Letters represent tide gauges in both panels (with red capital lettering denoting the tide gauges shown in Figs. 3, 5 and 8—A: brest (48.38° N, 4.49° W); B: Newlyn (50.10° N, 5.54° W); c: Malin Head (55.37° N, 7.33° W); d: Lerwick (60.15° N, 1.14° W); e: Aberdeen (57.14° N, 2.08° W); f: North Shields (55.01° N, 1.44° W); g: Immingham (53.63° N, 0.19° W); h: Lowestoft (52.47° N, 1.75° E); i: Zeebrugge (51.35° N, 3.20° E); J: Den Helder (52.96° N, 4.75° E); K: Delfzijl (53.33° N, 6.93° E); l: Esbjerg (55.46° N, 8.44° E); M: Hirtshals (57.60° N, 9.96° E); N: Smögen (58.35° N, 11.22° E); o: Bergen (60.40° N, 5.32° E). The red outline shows the extent of the region we included in regional mean statistics of the NWS

2.2.2 Satellite altimetry

We compare the CO6 simulated SSH mean and variability to satellite observations. We interpolate the altimetry products from their native resolution onto the CO6 model grid with bi-linear interpolation. As the satellite altimetry products are not able to resolve features with spatial wavelengths less than ~ 180 km (Legeais 2018), we smooth model and satellite product anomalies (after removing the domain mean) by convolving with a uniform filter of 13-by-13 grid boxes (~ 90 km).

The CO6 SSH mean (averaged from the last 200 year of the simulation) is compared to the AVISO Mean Dynamic Topography (MDT) product (0.25° horizontal resolution), estimated from the period 1993–2012 (Rio et al. 2014). MDT indicates the average strength of the geostrophic currents (Hermans et al. 2020a submitted). The CO6 interannual SSH variability is compared to that of the Copernicus Climate Change Service (C3S) Sea Level Anomaly (SLA) product (Legeais et al. 2018). We compute annual means (1993–2018) of the (daily, 0.25° horizontal resolution) C3S SLA product, linearly detrend, and calculate the temporal standard deviation. As the C3S SLA is 25 years long, we calculate the running temporal standard deviation of (linearly detrended) 25-year section of the last 200 years of model data. This gives a distribution of standard deviations, and we report the median, and 5th/95th percentiles of the distribution – these are compared to the C3S SLA Climate Change Initiative (CCI) (altimetry) data in Fig. 6.

2.3 Methods

We make use of the statistics outlined by Taylor (2001) to compare spatial patterns: the correlation (r) and the relative standard deviation (rsd). r tells us about the spatial similarity of the patterns. rsd compares the amplitude of the two patterns. Furthermore, following Taylor (2001), these two statistics can be used to calculate the (relative centred) RMS.

2.3.1 Analysis of modelled SSH

This study is focused on interannual sea-level variability, and so we analyse annual means from the model. Our standard analysis procedure is to extract the hourly mean values of SSH from the model, remove the tide with a Doodson filter (Pugh 1987), and calculate monthly and annual means from the residual. Given that the 360-day model calendar is made up of 12 30-day months, the annual mean can easily be calculated as an unweighted mean of the monthly means. Comparisons with (360-day) annual means of daily mean data show that the tidal aliasing is negligible on the annual timescale.

When comparing to tide gauge records, we take the nearest sea grid box to the tide gauge.

When we assess the fraction-of-variance of the SSH from the local steric and bottom pressure terms (Fig. 7), we use annual means calculated from the daily means, as the local steric and bottom pressure terms were not output at an hourly frequency.

2.3.2 Spectral analysis

We use spectral analysis to investigate SSH variance at different frequencies. The spectral analysis divides the records into sections (using a Hanning window), performs a discrete Fourier transform on each segment and averages the results. We also provide a crude estimate of the confidence of the modelled spectra, by providing the range of the individual spectra (as a grey shading around the modelled spectra (bold black line) in Fig. 5). We use a window of 64 years, with a 50% overlap (e.g., years 1–64 are the first segment, years 33–86 are the second segment, years 65–128 are the third segment, etc.). A 64-year window gives a frequency range of 64–2 year⁻¹. Increasing the amount of window overlap reduces the noise in the spectra, but also results in less well-defined spectral peaks. Given our relatively short timeseries, the 50% overlap provides a good compromise between accurately estimating signal power without over-counting any of the data. We use a boot-strapping technique adapted from Zwiers and von Storch (1995) to estimate significance of peaks and troughs within our spectra. To do this, we estimate the mean, variance and lag-1 correlation coefficient from the observed and model-simulated tide gauge. We then generate 1000 simulated AR(1) stochastic process random time-series with the same mean, variance and lag-1 correlation coefficient as the data. We calculate the spectra for each of these and find the 10th and 90th percentile for each frequency of the spectra. This provides significance thresholds for the spectral peaks and troughs (the dashed lines bounding the main spectra (bold) in Figs. 5 and 8)—any spectral trough or peak outside these thresholds is considered significant. The variance of the time series can be estimate by integrating below the power spectra, and this allows the variance to be separated into frequency bands. This is complicated in the method outlined above due to the sub-sampling and the Hanning window function. When we wish to estimate the variance at a particular frequency band (e.g., Fig. 10), we undertake a spectral analysis without subsampling, or using window functions (i.e., we use a periodogram).

Fig. 5

Power spectral density of the annual mean time-series for the observed (green) and model-simulated (black) tide gauge data. We use 64-years Hanning windows (with 50% overlap) for each Fourier transform, and so require time-series greater than 64 years to sample the lowest frequency (1/64 year⁻¹). When the observed time-series are less than this, we use a thinner green line to show the part of the spectra that is insufficiently sampled. The boot strapping technique gives significance limits as a pair of dashed grey/green lines—signifcant spectral peaks and troughs are outside these limits (see Sect. 2.3.2 for details). The grey shaded area around the model spectra (black) is a confidence interval (see Sect. 2.3.2 for details). The purple and olive shading show the spectral bands (8–20 years and 5–8 years respectively) chosen to capture the modelled spectral peaks and used in Fig. 10

2.3.3 Assessing the distribution of unforced sea level trends

Our long control simulation allows us to assess the distribution of unforced sea level trends, for a given period. We calculate the distributions with a running linear trend (an n-year trend calculated with data from n + 1 years). We convert this to an absolute change by multiplying by n (this allow different time periods to be easily compared). We use this to build a distribution for each grid box. As there are negligible SSH trends in the model, this distribution has a near zero mean, and we find they tend to be symmetrical (neutrally skewed). To show the maximum likely increase in sea level (for a given period) associated with natural climate variability, we produce maps of the 95th percentile of this distribution (Fig. 11). We note the 5th percentile map is very similar, but of the opposite sign. We also provide maps of the maximum likely trends in Additional Material Fig. 3.

2.3.4 Assessing how well a tide gauge represents the adjacent coast

In Sect. 4.4 we consider how representative tide gauge measurements are of the interannual variability of the adjacent coast. We extract the modelled coastline of four coasts: continental Europe (from the model boundary in Portugal to Denmark); the Scandinavian coast (from the Kattegat to the northern model boundary in Norway); Great Britain; and Ireland. For each grid box along these coasts we correlate the SSH time series with its neighbour and calculate the relative standard deviation of the two time-series. If the correlation and relative standard deviation meet some predefined threshold (i.e., r > 0.95, |rsd − 1.0|< 0.15), we move to the next neighbouring grid box, and continue until one of the thresholds is not met. This process is applied in both directions, to give the span of coast that can be represented by the tide gauge. We calculate the distance of this span by summing the grid box dimensions.

We can use this data to assess the distribution of values for each coast, and the span of coast represented by each tide gauge. As the shelf loading term (see Sect. 4.3) is the same everywhere, it increases the length-scales considerably. We therefore focus on the redistribution term, to give a conservative estimate. We tabulate the results presented in Fig. 12 in Table 12 and Table 13 in additional materials, and additionally give the equivalent values for the full SSH (loading and redistribution) term.

3 Model simulation evaluation

The main aim of this section is to assess the model’s ability to reproduce the observed characteristics of coastal sea-level variability, using the tide-gauge and satellite observations. However, in the additional materials, we include an evaluation of the fundamental behaviour of the shelf seas model (in terms of temperature and salinity) which gives further credibility to the modelling system. This is summarised in Table 3.

Table 3 Summary of model evaluation

3.1 Evaluation of tide gauge interannual variability

The focus of our model evaluation is the comparison with interannual tide gauge variability. There are several processes that are either not resolved, or not represented in the model that influence the observational timeseries. These include harbour seiches, geophysical processes (such as tectonic motions), glacial isostatic adjustment and land subsidence (e.g., Douglas 2001). In addition, tide gauges can be subject to instrument errors and biases (e.g., changing datums). These possible errors must be taken into consideration when comparing observed tide gauges to model output.

The first stage is to compare the 32-year high pass filtered annual mean time-series from the observed and model-simulated tide gauges. These time-series are plotted in Fig. 3 with their normalised distributions, and sample standard deviation. The observed tide gauges tend to have a greater interannual variability than the models, and a greater observed kurtosis (fatter tails in the distributions). Most observed and modelled distributions tend to be symmetrical (with absolute skewness < 0.5) and the median (denoted by a bold horizontal green/black line against the distributions in Fig. 3) near the mean (zero). Tide gauges that do have skewed distributions often have large gaps within the tide gauge record (e.g., Esbjerg) or low frequency variability (e.g., Smögen)—this could reflect missing processes within the model. The observations are often platykurtic (kurtosis > 0 having fatter tails), while the modelled tide gauges tend to be mesokurtic to leptokurtic (kurtosis < 0). The tide gauge with the shorter observed records tend to have greater kurtosis—this may reflect a few outlier values having a disproportionate affect, e.g., Malin Head.

There is generally a good qualitative agreement in the shape of the modelled and observed distributions, particularly for the longer tide gauge records. For the shorter records the irregular distribution reflects paucity of data (e.g., Lerwick, North Shields), or possible inconsistencies in the data record (e.g., Lowestoft). For the longer tide gauge records there is an excellent qualitative agreement between observed and modelled distribution at Newlyn, Brest and Den Helder. The observed tide gauged at Delfzijl and Hirtshals appears to have secondary peaks, which are not captured by the model. This could reflect missing modelled processes or issues with the tide gauge data.

The second stage of evaluation against tide gauges is to consider the spatial pattern of the interannual-to-decadal variability. We compare the sample standard deviation of the observed and model-simulated 32-year high pass filtered tide gauge in Fig. 4. We find a good spatial agreement between the two datasets, with the model able to reproduce the pattern of shelf sea-level variance. The two spatial patterns of the observed and model-simulated tide gauge variances are significantly correlated (spatial Pearson’s correlation coefficient r = 0.78). The magnitude of the pattern (its standard deviation) is greater in the observed tide gauges (with a rsd = 0.59), reflecting the model’s underestimation of the tide gauge variability by ~ 33% (the linear slope coefficient m = 1.333). This is expected due to limitations in the represented processes and around the deep low pressures and peak wind speeds. Our map of tide gauge interannual variances agrees visually with Wahl et al. (2013).

3.2 Evaluating against tide gauges power spectra

We have shown that the model is able to simulate the magnitude and spatial pattern of the observed tide gauge interannual variability and the distribution and higher statistical moments at selected locations. We now assess the model’s ability to simulate the spectra of this variability at these selected locations. Details of the years used in the spectral analysis is given in Table 2.

We find similar peaks in our observed spectra to Unal and Ghil (1995) at Brest (they find peaks at 32, 17.1 and 12.2 years), Newlyn (10.7 years), Smögen (18.03 and 10.2 years) and Hirtshals (32 year).

Overall, there is a good agreement between the model-simulated and observed tide gauge spectra with many shared characteristics and features. All show a decrease in energy towards higher frequencies (note the logarithmic scale). The observed tide gauges are typically within the (crude) uncertainty estimates of the model-simulated tide gauges. Most locations show an energy trough at about the 8-year period (sometimes towards 10 years) in both the observed and modelled tide gauge (with the exception of observed tide gauge record at Smögen and Brest), which is generally significant (with the exception of the modelled Newlyn tide gauge record). The observed spectra at Brest and Smögen have a similar trough separating two peaks, but this is at a lower frequency (~ 14 years)—the modelled spectra have the trough at ~ 9 years. We have less confidence in the estimates for the lowest frequencies, as less cycles can be captured. The observed and modelled ~ 8-year trough (~ 14-year for the observed spectra at Brest and Smögen) typically separates two peaks. The higher frequency peak is significant for both observed and modelled tide gauges for most locations. For the modelled spectra this peak typically extends from 5 to 8 years, whereas the observed spectral peaks are often extending to higher frequencies.

At frequency lower than 8 years, there is typically a broad, insignificant peak. For Newlyn, Den Helder and Delfzijl, this is a defined peak in the observations, which is significant for Newlyn and Den Helder. Our observed records at Hirtshals and Smögen (57 and 55 years respectively, see Table 2) do not capture two cycles of a 32 year period, so we are less confident of behaviour of the spectra beyond this apparent peak. Likewise, Brest has a significant observed peak which begins at 32 years (in agreement with Unal and Ghil (1995))—it appears to extend to lower frequencies but as our 64-year frequency estimate is only sampled once in the 83 year record we have less confidence in this value. We note that the Brest tide gauge observational uncertainty reported by Wöppelmann et al. (2008) does not apply here as we use Brest data from 1861–1943. Like Hirtshals, Smögen has a lower frequency peak, but due to the observed record length, it is unbounded at lower frequencies.

The modelled spectra all have broad, insignificant low frequency peaks bounded by this ~ 8-year trough. Brest and Newlyn have similar peaks between 8 years and 32 years, and a suggestion of a peak at lower frequencies. The North Sea (Den Helder and Delfzijl) and Skagerrak (Hirtshals and Smögen) tide gauges have similar low frequency peaks, that extend from a (near significant) trough at 8–64 years, with local maxima at 12 year and ~ 20 years.

In Sect. 4.2 we separate our data into a low and high frequency component. We use the two spectral peaks and the ~ 8-year trough to define two spectral bands: 5–8 years, and 8–20 years. These are highlighted (with shading) in the subsequent spectra plot (Figs. 5 and 8).

3.3 Evaluation of SSH against satellite altimetry data

We compare the Ctrl mean SSH to the AVISO MDT in Fig. 6a, d. Both the model and altimeter have very similar spatial patterns (with significant spatial correlations = 0.95, rsd = 0.98 across the domain), and both fields show a higher mean SSH in the German Bight, Skagerrak/Kattegat, in some coastal regions, and to the southwest of the NWS (in the open ocean). All show lower mean SSH to the north of the NWS (i.e., the Norwegian Sea, north of the Greenland-Scotland ridge), although this region is larger and more intense in the AVISO MDT (Fig. 6a). The MDT reflects the mean geostrophic currents. On the NWS, the Ctrl captures the mean circulation well (e.g., Turrell et al. 1992; OSPAR 2000, not shown), and this is reflected in the modelled mean SSH field. For example, the zero contour on Fig. 6a runs parallel to the Scottish coast, before entering the North Sea by the Orkney Islands. This represents this main North Sea inflow. This contour then crosses the North Sea as the Dooley current, before turning south into the Norwegian Trench, and returning as the Norwegian Coastal Current. Likewise, the 100 mm contour in the southern North Sea represents the secondary, English Channel North Sea inflow. This all matches with the established NWS circulation, giving further confidence in Ctrl. The AVISO MDT (Fig. 6d) also matches this circulation, although doesn’t do quite so well in places (e.g., the north-eastern North Sea).

Fig. 6

Comparison of model-simulated (Ctrl) and satellite observed MDT and SSH variability. Left column: a model-simulated MDT anomaly and d AVISO MDT anomaly (with the spatial mean removed from both, mm). Centre and right columns: SSH variability (standard deviation, mm). Ctrl is broken into 25-year samples (with 50% overlapping), to build a distribution of interannual standard deviation of a 25-year period—the range of this distribution is presented as the 5th and 95th percentile (e, b respectively) about the median (c). The C3S Sea level Anomaly CCI (f) can be compared to the median (c) and can be expected to fall within distribution (5th–95th percentile range, e, b respectively)—where it does not (i.e., f is outside the range e, b), the C3S (f) data is stippled with red. Spatial correlation coefficients and relative standard deviation are given for the relevant observations (panel a and f). All fields are smoothed with an ~ 90 km length-scale. The inverse barometer effect is removed from the modelled data for this plot (it is included in all other analysis and visualisation within this study)

We compare the distribution of modelled SSH variability to the C3S SLA variability. There are many features of the spatial pattern of variability common to both datasets, with a NWS spatial correlation of r = 0.89 (rsd = 1.45) between the C3S and the median of the modelled distribution (Fig. 6c, f). In the open ocean, in both the observations and the model, there is enhanced sea-level variability in deeper water (e.g., the Icelandic basin, Rockall trough) and there is suppressed variability in the vicinity of the shelf break, particularly to the south west, i.e., bounding the Celtic Sea and the Armorican shelf. On the shelf the greatest sea-level variability is in the German Bight. The spatial pattern of the C3S SLA variability has a larger amplitude than Ctrl. On the NWS, the C3S SLA observed variability is largely within the modelled distribution of variability. In the open ocean, the observed variability is generally larger than the distribution of modelled variability.

We have performed the same analysis on the GC3.0 SSH fields (Additional Material Fig. 2) to show the impact of downscaling. The main reason for the downscaling is to include physical processes important to the NWS SSH that are absent in the GC3.0, and so assessing the impact of this downscaling is not an aim of this study. GC3 approximately captures the general circulation of the NWS, but with some important differences. For example, the main North Sea inflow, and northern North Sea circulation is not as well represented, nor is the English Channel North Sea inflow. The C3S SLA interannual SSH variability is within the CO6 modelled distribution for most of the NWS (Fig. 6f), whereas the is it largely outside the GC3.0 distribution on the NWS (Additional Material Fig. 2f).

We have undertaken extensive model evaluation and have shown that the model is suitable to simulate near present day sea-level interannual-to-decadal variability around the UK and within NWS. We now use the model to investigate some of the sources of the NWS sea-level variability and consider its implications for sea-level projections.

4 Interannual SSH variability

4.1 Relative importance of local steric and bottom pressure terms

Firstly, we investigate how much of the local steric and non-steric (bottom pressure or manometric) components contribute to the total variance (see Sect. 11.2 in the methodology appendix for details of how they are diagnosed from the model). The bottom pressure component is related to water moving on and off the shelf (shelf loading) and moving around within the NWS (redistribution). The local steric effect relates to the local expansion of the water column by warming and freshening. As this expansion is a depth-integrated effect, its possible magnitude decreases with the water depth, and so is relatively small over most of the NWS. Steric anomalies in the open ocean can propagate onto the shelf as mass signals, and so are reflected in bottom pressure (e.g., Landerer et al. 2007; Bingham and Hughes 2012). Within our study we are not able to distinguish between the remote steric effect and the barotropic sea level.

Following Roberts et al. (2016), we calculate the Fraction-Of-Variance (FOV) to show how much of the Ctrl SSH variance is explained by changes in local bottom pressure (i.e., column-integrated mass) and local steric changes (i.e., column-integrated density) (Fig. 7, see Sect. 11.3). On the shelf, the variance associated with bottom pressure dominates the SSH variance (Fig. 7), accounting for 95% (when averaged over the shelf). Off the shelf, particularly to the west, the bottom pressure term is smaller. There is compensation between the steric and barotropic components where the SSH variance is smaller than the sum of the variance of the two components—this corroborates Roberts et al. (2016). On the shelf, the local steric change accounts for 13% of the total variance (when averaged over the shelf). Off the shelf, the local steric component becomes much more important. Since our primary focus is on-shelf SSH variability, and modulation of coastal flood risk, we disregard the local steric effects for the remainder of this study. We note that the FOV varies with season (not shown). In summer (June–August), the steric (bottom pressure) accounts for 14% (91%) of the shelf mean SSH variance, while the other seasons it is closer to 6% (97%).

Fig. 7

Fraction of variance (FOV) of the interannual sea-level variability associated with a the bottom pressure and b the steric effect. Negative values (grey) suggest where there is compensation between the two components

4.2 Isolating sources of interannual SSH variability associated with the model boundary conditions

We can isolate the NWS interannual SSH variability associated with lateral ocean and surface atmospheric boundary variability with a pair of sensitivity experiments: CtrlOcV and CtrlAtV (Table 4, holding the atmosphere to climatology (with a seasonal cycle) in CtrlOcV and ocean to climatology in CtrlAtV). As with Ctrl, the river and Baltic exchange boundary conditions are forced by a seasonal climatology in CtrlOcV and CtrlAtV. We find that the interannual SSH variance of the two sensitivity experiments combine linearly, with a shelf mean error and covariance term that combine to less than 0.1% (see spatial means of SSH spatial standard deviations in Fig. 9a, e, and f). We analyse CtrlOcV and CtrlAtV spectrally (c.f. Fig. 5), to assess how these drivers affect different parts of the spectrum (Fig. 8). We then assess the spatial patterns of variability of these different drivers in Fig. 9.

Table 4 Descriptions of model runs

Fig. 8

Comparison of power spectral density of the sensitivity experiments. See Fig. 5 for details. For each (model-simulated) tide gauge the black line represents Ctrl (the same line as in Fig. 5). The blue (red) line represents CtrlOcV (CtrlAtV), where the atmosphere (lateral ocean) forcings are kept to a climatology, and so all the forced variability comes from the ocean (atmosphere) boundary conditions

Fig. 9

The left hand column is the total variance (var_ssh: a, e, i), the second and third column is the variance associated with the resdistribution (var_red: b, f, j) and loading term (var_lod: c, g, k) respectively, with (two times) their covariance in the right hand column (2 × cov_red,lod: d, h, l). The total variance (var_ssh) is the exact sum of the variance associated with the resdistribution and loading and two times their covariance (var_red, var_lod and 2 × cov_red,lod). The variance from Ctrl (upper row) is the approximate sum of the variance CtrlOcV (middle row) and CtrlAtV (bottom row), with any difference being due to non-linear interactions. White contours match the values of the colourbar ticks, and the bold white contour in the covariance panels represents 0 mm². The NWS mean (and 5th–95th percentile range) is given in Table 5 for each panel

There is a large decrease in energy towards high frequency (> 1 year⁻¹) when the atmosphere is held to a climatology (CtrlOcV, blue in Fig. 8). This suggests that little of the high-frequency SSH variability on the shelf comes from the ocean boundary conditions. This is corroborated by the convergence of Ctrl and CtrlAtV (the black and red lines) at frequencies higher than about ~ 8 years (~ 5 years for Newlyn and Brest). At frequencies lower that this, the red and black lines are parallel, but the SSH variability associated with atmospheric variability is notably lower. Conversely, at very low-frequencies (< 0.02 year⁻¹, ~ 50 years), there is a suggestion that Ctrl converges with CtrlOcV (particularly in the German Bight sites, Fig. 8c–f), implying that at the lowest frequencies, the ocean boundary conditions are the dominant source of NWS SSH variability. The ocean driven climate variability such as the Atlantic Multidecadal Variability (AMV) and AMOC could lead to SSH variability at these low frequencies (Häkkinen et al. 2013). These results agree with Dangendorf et al. (2014, 2015) who separated the tide gauge variability associated with local atmospheric forcings, and showed they dominated the higher frequencies, while the remaining residual term was important at lower frequencies.

The variability associated with the atmosphere (CtrlAtV, red) tends to have similarly located peaks (~ 13–16 and ~ 7 years) and troughs (~ 8 years) to Ctrl (described in Fig. 5), although the lower frequency peak (~ 13–16 years) has less power. The spectral peak associated with the ocean variability (CtrlOcV, blue) is very different from Ctrl (black). There is a very defined broad (significant) peak at 7 years between two (significant) troughs at 16 years and ~ 5 years. The spectra of CtrlOcV are similar between tide gauges (not shown), suggesting very little spatial variability and suggesting (but not proving) signal coherence.

The spatial patterns of the SSH variance (Fig. 9a, e, i) give further insight into Fig. 8. When looking at CtrlOcV SSH variability (Fig. 9e), we see that there is negligible spatial variance on the shelf (reflecting the similarity of the spectra between tide gauge locations). We also investigate the coherence of the CtrlOcV NWS SSH signal, by looking at the temporal correlations between grid boxes (not shown) and find high correlations (typically r > 0.95), and near unity relative standard deviations (typically 0.95 < rsd < 1.05). This indicates that the CtrlOcV SSH changes are coherent, with the whole sea level rising and falling as a level surface across the NWS, on the interannual timescale.

CtrlAtV NWS SSH variance (Fig. 9i) has a much stronger spatial pattern. The greatest variance is in the eastern North Sea (in the German Bight and along the Danish coast), and around Scotland, in agreement with Hermans et al. (2020a submitted). We find that NWS SSH variance associated with the oceanic and atmospheric boundary variance combines almost linearly in Ctrl (Ctrl Fig. 9a, e, i). As there is little spatial pattern in CtrlOcV, the spatial patterns of Ctrl and CtrlAtV NWS SSH variance are highly correlated (r = 0.98, rsd = 1.04).

4.3 SSH variability associated with loading and redistribution

NWS SSH variance can be separated into components. The total amount of water on the shelf can vary—we refer to this term as “shelf-loading”. The redistribution of this water within the NWS can also lead to local SSH variance. We separate SSH into a shelf-loading time-series (ssh_lod, a proxy for a time-series of the total volume on the shelf, which is spatially homogenous) and a volume redistribution term (ssh_red; which, by construction, always spatially averages to zero). To separate these terms, we first remove the modelled temporal mean SSH value from each grid box to give the SSH temporal anomaly (ssh_anom):

$${ssh}_{anom(x,y,t)}={ssh}_{\left(x,y,t\right)}-\frac{\sum_{t=0}^{T}{ssh}_{\left(x,y,t\right)}}{{n}_{t}}$$

(1)

We then calculate a time-series of the ssh_anom spatial mean which is the SSH change associated with the change in the total water mass on the NWS. We call this the shelf loading term (ssh_lod).

$$ss{h}_{lod\left(t\right)}= \frac{{\sum }_{x}{\sum }_{y}{ssh}_{anom(x,y,t)}}{{n}_{x}{n}_{y}}$$

(2)

The redistribution term is the difference between the anomaly and the shelf loading term:

$$ss{h}_{red\left(x,y,t\right)}={ssh}_{anom(x,y,t)}-ss{h}_{lod\left(t\right)}$$

(3)

The variance of these two terms, with their covariance, combine (exactly) to give the variance of the original time-series.

$${{\sigma }^{2}}_{{ssh}_{anom(x,y,t)}}= {{\sigma }^{2}}_{ss{h}_{red\left(x,y,t\right)}}+{{\sigma }^{2}}_{ss{h}_{lod\left(t\right)}}+{2\times cov}_{\left(ss{h}_{red\left(x,y,t\right)},ss{h}_{lod\left(t\right)}\right)}$$

(4)

We refer to these values as:

$${var}_{ssh}={var}_{red}+{var}_{lod}+{2\times cov}_{red,lod}$$

(5)

As the ssh_red is the anomaly term around ssh_lod (3) and so has a zero NWS spatial mean, the covariance term also has a spatial mean of zero across the shelf.

When comparing the power spectra of the shelf loading timeseries for CtrlOcV with the modelled SSH spectra at any of the tide gauges (not shown) we find an excellent agreement with the dominant 5–16-year band well represented. The spectra of the Ctrl also generally agree, capturing both the 5–8- and 8–20-year bands.

The spatial patterns of the loading and redistribution variances and their covariance for the three runs are shown in Fig. 9, with the NWS mean and 5th and 95th percentile values given in Table 5. The variance of the loading terms is, by design, spatially homogenous (Fig. 9c, g, k), and is associated with both the oceanic and atmospheric variance (Fig. 9g, k), although the CtrlAtV is greater. These two terms add up almost linearly in Ctrl, suggesting the non-linear interaction is small. Volume redistribution on the shelf is almost exclusively driven by the atmospheric variability (c.f. Fig. 9f and j). While the variance of the redistribution terms is smaller than that of the loading term (~ 20% of the shelf mean variance in Ctrl, ~ 35% in CtrlAtV), it has a substantial spatial pattern (Fig. 9j). The greatest redistribution variance occurs in the eastern North Sea (in the German Bight and along the western coast of Germany and Denmark), while there are slightly raised values in the southern Celtic Sea.

Table 5 Summary statistics (the NWS mean, and 5th and 95th percentile values) associated with Fig. 9

The covariance of the shelf-loading and redistribution terms shows where these terms reinforce or interfere with one another. It is driven by atmospheric variability (as the CtrlOcV redistribution term is negligible, its covariance is also negligible), and (by construction) has a spatial pattern with a regional mean of zero. The covariance term (Fig. 9d, l) is greatest in the eastern North Sea (perhaps centred further north than in the redistribution term) and around Scotland, where it is the dominant source of variability. It is negative in the Celtic Sea, English Channel and the south-west North Sea. This term is particularly important in two regions: the Celtic Sea and around Scotland. The covariance tends to cancel out the effect of the redistribution variance in the Celtic Sea. Around Scotland, the total variance (var_ssh) is high, while the redistribution variance is small, and so the signal comes from the covariance.

To further investigate these spatial features, we return to the spectral bands identified in Fig. 8. Since the loading terms and CtrlOcV have little spatial pattern, we focus on Ctrl and CtrlAtV, the total variance (var_ssh), the redistribution variance and the covariance term (although provide summary statistics for all terms in Table 6). We separate these terms into 5–8 years (7-year peak) and 8–20 years (13-year peak) frequency bands (Fig. 10). We effectively split the variance from each panel in Fig. 9 into the two frequency bands (as described Sect. 2.3.2, with the covariance term calculated by subtraction). We remove the second row and third column from Fig. 9 as they have minimal spatial patterns. The remaining panels (Fig. 10) do not quite add to give Fig. 9, as we have removed the variance associated with frequencies higher than 5 years, and lower than 20 years. The Ctrl total variance (var_ssh) is separated fairly evenly between the two bands (with area means of 127 mm² and 131 mm² for 5–8- and 8–20-year frequency bands). While this having a near unity ratio is relatively arbitrary (it reflects the size of the peaks in Fig. 8, and also the width of the spectral bands) it allows us to compare the other terms. The CtrlAtV is much lower at the lower frequency (99 mm² and 68 mm² respectively), reflecting the CtrlAtV lower energy peak begin weaker than Ctrl in Fig. 8. We noted that the CtrlAtV 5–8-year peak was similar to Ctrl and for the North Sea tide gauges in Fig. 8 but was weaker in Brest and Newlyn. This is clear when comparing Fig. 10a and b, with Ctrl having much higher energy to the south and west of Great Britain. Conversely, Table 6 shows that CtrlOcV is has more energy at the lower frequency band (36 mm² and 44 mm²).

Table 6 Summary statistics (the NWS mean, and 5th and 95th percentile values) associated with Fig. 10

Fig. 10

Frequency dependence of spatial patterns of the total variance (var_ssh), the redistribution variance (var_red) and the covariance term (cov_red,lod). Panels from Fig. 9 (a, b, d, i, j and l) with substantial spatial pattern (excluding the loading term (var_lod), and CtrlOcV), are separated into a high frequency (5–8 years) and low frequency (8–20 year) band (these bands are shown in the spectra in Fig. 5). The rows give the total variance (var_ssh: a–d), the variance associated with the redistribution term (var_red: e–h) and (2 times) the covariance term (cov_red,lod: i–l). First and second columns are the high frequency (5–8 year) component (a, b, e, f, i, j), while the third and fourth give the low frequency (5–20 year) component (c, d, g, h, k, l). The first and third column are from Ctrl (a, c, e, g, i, k), while the second and fourth columns are CtrlAtV (b, d, f, h, j, l). Summary statistics (NWS mean, 5th and 95th percentile values) are provided in Table 6. Note that the colour ranges are different for each pair of panels

Ctrl and CtrlAtV have very similar patterns (and magnitudes) of redistribution variance at both frequencies (reflecting spatial coherence of CtrlOcV). The patterns of the variance associated with the redistribution term (for both Ctrl and CtrlAtV) are similar at the two bands, but with some subtle differences. The low frequency redistribution variance (Fig. 10g, h) is greater (than the high frequency component) around northwest Scotland, in the southern North Sea along the pathway of the south North Sea Water (OSPAR 2000; which flows north-westward from East Anglia, towards the Skagerrak) and in the vicinity of the seasonal tidal mixing front (Brown et al. 1999).

When looking at the covariance term, the difference between frequency bands is clearer. Although the general pattern is similar, there are important differences. The covariance in the eastern North Sea is much greater at higher frequencies (Fig. 10i, j compared to Fig. 10k, l), and the latitude of the peak values also changes. At higher frequencies (Fig. 10i, j), the covariance is greater along the east coast of the Jutland Peninsula, whereas at lower frequencies (Fig. 10k, l) there is a clear local maxima at the northern end of the Jutland Peninsula. This is collocated with a reduced low frequency redistribution variance, suggesting a change in the correlation rather than a change in the variance. The total variance around Scotland is slightly greater at lower frequency (Fig. 10a, c) whereas the CtrlAtV covariance is slightly greater at higher frequencies (Fig. 10j, l). Figure 10 also illustrates how the covariance and redistribution terms cancel out in the Celtic Sea and western English Channel, and this is predominantly a high frequency atmospheric process (Fig. 10f, j).

4.4 Internal variability in coastal sea level: trends and coastal coherence

Projected twenty-first century sea-level change is often illustrated as a set of smooth curves (e.g., Church et al. 2013). Over a period of several years, trends associated with temporal variability can dominate over the projected long-term trends. Palmer et al. (2018) showed that over the next decade sea-level uncertainty is dominated by internal variability rather than model-structural uncertainty and emission uncertainty (Fig. 1)—users interested in sea level over the next few decades must consider sea-level variability. Traditionally, long tide gauge records have provided the best estimates of sea-level variability on a range of timescales. These estimates are spatially sparse compared to the scale of coastal engineering and development projects, and it is unclear how representative each estimate is for the surrounding coastline. To aid interpretation of the tide gauge variability estimates, we analyse the spatial range of this SSH variability.

We assess the maximum likely absolute sea-level change associated with natural variability for 4 different periods. We have identified two peaks at 7 and 13 years in Fig. 5 (associated with the spectral bands used in e.g., Fig. 10). When looking at the spectra of all the tide gauges in Fig. 4 (not shown) we find additional peaks common to many tide gauges at ~ 4 years and ~ 21 years. The exact values are in part related to the resolution of the spectral analysis. Following the methodology outlined in Sect. 2.3.3 we show the greatest likely (95th percentile) change in SSH in Fig. 11 (with summary statistics included in Table 7). To put these into the context of the present-day rate of sea level change (3.2 mm/year Church et al. 2013), and projected end-of-century rate (12 mm/year, RCP8.5, Church et al. 2013), we give the equivalent change in sea level at 3 and 12 mm/year for each period in the respective panels in Fig. 11. We also include an equivalent figure of the maximum likely trends in the Additional Materials (Additional Material Fig. 3), while the summary statistics are also included in Table 7. When averaging over the shelf, there is a possible 58 mm increase in regional sea level over a 4-year period (58.9 mm when averaged around the coastline of Great Britain). The greatest values are in the eastern North Sea for all time periods apart from 13 years, where the values around Scotland have a similar magnitude. The shelf-mean of the trends on Fig. 11 decrease with time period.

Fig. 11

Maximum likely (at the 95th percentile) absolute increase (trend times number of years) in unforced sea level over 4, 7, 13 and 21 years. The NWS mean (outlined by the red contour) is given in each panel, and the coastal mean for the coast of Great Britain, Ireland, Europe (from the model domain boundary in Portugal, to the Skagerrak), and Scandinavia (from the Skagerrak to the model domain boundary in Norway) is given in Table 7. The maximum likely trends are given in Additional Material Fig. 3

Table 7 Mean value of the 95th percentile unforced sea level change (and trend), along different coasts, for each period. From Fig. 11

There is a large degree of spatial coherence in the sea-level variability, suggesting that the simulated tide gauge variability estimates are generally representative of a much greater length of coast. We follow the methodology outlined in Sect. 2.3.4 to assess the length of coast that each tide gauge (and each coastal grid box) represents. We use a correlation threshold of r = 0.95 and expect the standard deviations to be within 15% of one another (Fig. 12). As the shelf loading term is spatially invariant, its inclusion increases this length substantially, and so here we focus on the redistribution term, but provide both values in Table 12 and Table 13 in additional materials. We note that, when using the full SSH time series (ssh_red + ssh_lod), the length scales are increased by low pass filtering the data, but this is much less important when looking at the redistribution SSH time series.

Fig. 12

The natural variability of each tide gauge (and coastal grid box) is correlated with the adjacent grid boxes along the coast. Here we show the length of coast that can be represented by the interannual variability of a given tide gauge (with r > 0.95, |rsd − 1|< 0.15) as coloured sections on the map. These sections are given some transparency, to show where neighbouring sections overlap. These coast lengths are tabulated in Table 12. For each of the modelled major coasts (those of Europe, Scandinavia, Great Britain and Ireland), we give the distributions of all coastal grid boxes (with 5th, 50th (median) and 95th percentiles, and mean (bold) denoted with vertical lines). The summary statistics of these distributions are tabulated in Table 13 in additional materials

These values are much higher when using the full SSH timeseries (ssh_red + ssh_lod)—for example the median value for the coast of Great Britain increases to 1177 km and Europe to 905 km. In the real world, other processes not represented in our model (such as harbour seiches, tectonic movement, lard subsidence etc.) may reduce this coherence.

5 Discussion

We outline a series of aims in the introduction. Our preliminary aim is to assess the model’s ability to simulate interannual variability, and then to describe the behaviour of the interannual variability, and finally to provide some end-user guidance. Here we discuss the results in the context of these aims.

We assess the model simulated sea level, temperature and salinity. The temperature and salinity are within the range of other similar downscaling studies. The SST biases are 0.4 °C when averaged over the NWS (up to 0.9 °C in the spring/summer North Sea), with the OSTIA observations being within two standard deviations (the modelled interannual variability) of the model over most of the NWS (apart from the summer eastern North Sea). The upper salinity biases are typically < 0.5 psu. As the focus of the study is sea level, we compare our simulations against tide gauge records and satellite altimetry. There are several important processes not included in our model simulation which will impact the tide gauge observations. This puts an upper limit on possible observation model agreement. We show that there is generally a good agreement between the modelled and observed sea level distributions (Fig. 3), and interannual variability when compared spatially across the NWS (spatial correlation of r = 0.78, Fig. 4). We select six tide gauges with particularly long records to compare the spectra of observed and modelled sea level variability. There are many common features between the observed and modelled tide gauges, including a trough at about 8 years separating two peaks of similar spectral width and power. We show that the modelled MDT (indicating geostrophic currents) has an excellent agreement on the NWS with the altimetry data. We also find downscaling improves the qualitative agreement between modelled MDT and NWS circulation (when compared to the GC3.0 MDT). The altimetry observed interannual variability is within the modelled distribution over most of the NWS, although is greater in the open ocean. From this extensive model evaluation, we conclude that the model is able to provide useful simulation for the analysis of SSH interannual variability on the NWS.

The local steric term is small compared to the bottom pressure term, as expected given the depth of the NWS, corroborating other studies (e.g., Landerer et al. 2007; Roberts et al. 2016). When the ocean adjacent to the NWS warms and expands (steric sea-level rise), it induces a horizontal pressure gradient, which is balanced by water flowing onto the NWS (e.g., Landerer et al. 2007; Yin et al. 2010; Bingham and Hughes 2012; Richter et al. 2013)—this occurs in non-Boussinesq and Boussinesq models alike (Griffies and Greatbatch 2012). The remote steric effect is therefore incorporated in the bottom pressure records in our simulations, and can be an important contribution, especially on longer (decadal) timescales (Dangendorf et al. 2014). Chen et al. (2014) were able to isolate this term, by running additional barotropic runs—we did not do this. The NWS has an important baroclinic circulation component (e.g., Hill et al. 2008), so this approach may introduce other errors.

We use sensitivity tests to isolate the variability associated with the atmospheric and oceanic lateral boundary conditions (referred to as “CtrlAtV” and “CtrlOcV” respectively). The atmospheric and oceanic terms combine almost linearly, with non-linear interaction terms < 0.1% when averaged over the shelf. The sea-level variability on the shelf caused by the ocean boundaries has little spatial pattern and is coherent implying the sea level across the shelf rises and falls as one, as suggested by the results of Dangendorf et al. (2014). This changes the total amount of water on the NWS, without changing its spatial distribution—we term this “shelf loading”.

The sea level of the North Sea and the wider NWS is particularly affected by changes in the atmosphere (Dangendorf et al. 2014). We find that the atmosphere driven variance is typically greater than the oceanic driven variance. We do not further separate the oceanic boundary forcings into the steric and barotropic terms, or the atmospheric boundaries into wind, pressure and heat/moisture fluxes (e.g., Hieronymus et al. 2017; Hermans et al. 2020a submitted). Hermans et al. (2020a submitted) ran similar sensitivities tests to this study, with additional tests to show the variance associated with wind, pressure, and buoyancy fluxes. Their SSH variability associated with atmospheric variability had a similar spatial pattern to ours, and they showed that this pattern mainly comes from the wind forcing—the variability associated with the inverse barometer effect, and buoyancy fluxes was fairly homogeneous over the shelf. The inverse barometer effect was most important where wind was less dominant, and while the buoyancy fluxes had a smaller impact they should not be neglected—the balance between these terms depended on location. Hermans et al. (2020a submitted) also showed that the southern and western ocean boundaries were more important that the northern (and Baltic) boundaries.

The CtrlOcV SSH variance is spatially uniform and coherent. This is consistent with a sea level anomaly in the ocean adjacent to the NWS driving a fast barotropic wave adjustment on the shelf. Such an adjustment would be fast compared to the rate of sea-level change, hence would appear as a uniform change. The shelf slope bathymetry and current system limits the advective exchange between the shelf and the open ocean (e.g., Wakelin et al. 2009). Understanding the shelf sea level response to the wave and advective processes and their relationship to future changes to the shelf slope current is an interesting topic for future research. Changes in shelf sea level, particularly in the coastal regions, can affect the tidal propagation and can have an important impact on coastal sea level (e.g., Howard et al. 2019; Haigh et al. 2020).

The atmosphere variability leads to a loading of water on the NWS, but also redistribution of water within the NWS. In addition to CtrlAtV having 92% greater shelf loading than CtrlOcV, it also has an important redistribution term, which has a considerable spatial pattern. The separation of CtrlAtV (and Ctrl) into loading and redistribution leads to an important covariance term. The redistribution and covariance associated with CtrlAtV are the only terms that provide a spatial pattern in interannual sea-level variance on the NWS.

CtrlAtV covariance does not contribute to the (spatially averaged) shelf mean SSH variance (it averages to zero across the NWS) but is an important term in some locations. For example, the high sea-level variance around Scotland is due to the covariance, as the CtrlAtV redistribution is relatively small here. This term reflects the interaction of the water moving onto the shelf and water moving around the shelf, so is positive in places where the amount of water on the shelf increases (the loading term) when the local sea level is also rising relative to the shelf mean sea level (the redistribution term). It is therefore not surprising that the covariance term is positive around Scotland, as this is one of the main routes of the Atlantic water into the North Sea (Turrell et al. 1996). There is a large region where covariance < 0 in the Celtic Sea, suggesting the loading and redistribution terms are out of phase in this region. This suggests that when the sea level in the northern North Sea (and German Bight) is rising (relative to the shelf mean, loading term), it is falling in the Celtic Sea (relative to the shelf mean). Wind patterns that drive the loading term may have important regional structure linking these local relative sea level anomalies (relative to the to the shelf mean SSH i.e., the redistribution term). Exploring the linkages between the wind drivers, shelf circulation and NWS sea level response will be the subject to further research.

SSH changes arise from a number of different processes and mechanisms. When analysing the spectra from the observed and model-simulated tide gauges, we find distinct spectral peaks, some of which appear to vary between the different tide gauges, and others that appear spatially homogeneous (Fig. 8). We speculate that a spatial frequency decomposition may show different mechanisms working at different location and frequencies. We divide our variance terms into high and low frequency components (5–8 and 8–20 years), as shown in Fig. 10. For each frequency, we expect Ctrl to be the sum of CtrlAtV and CtrlOcV. The total variance of Ctrl has a similar pattern for both frequencies. The total variance of CtrlAtV is (relatively) much weaker in the low frequency band, reflecting the spectral characteristics of the tide gauges in Fig. 8. We therefore expect CtrlOcV to be greater in the low frequency band. The CtrlOcV variance is ~ 20% greater in the low frequency band (i.e., CtrlOcV var_ssh (= var_lod) is 36 mm² and 44 mm² in the 5–8 and 8–20 year frequency bands—Table 6), however this is not very clear in Fig. 8. The Ctrl and CtrlAtV redistribution terms are similar in magnitude and pattern at both frequency bands (Fig. 10e–h), therefore the spatial difference is largely due to the covariance term. The spatial pattern of the covariance term is frequency dependent in places. For example, the covariance term shows frequency dependence along the Jutland Peninsula. The pattern of the Ctrl var_ssh local maxima is similar at both frequency bands (Fig. 10a, c) in this region. At high frequency, the pattern of the redistribution and covariance terms are similar (Fig. 10e, i), however they are different at low frequency (Fig. 10g, k). At low frequency, the redistribution term is dominant to the south (Fig. 10g), while the covariance term is dominant to the north (Fig. 10k), reflecting the changing correlations. This is perhaps an example of different processes working at different frequencies. As well as being a region with large SSH variability, it is also the main pathway of water flowing into the Baltic Sea, and an important pathway into the Norwegian Coastal Current via the Skagerrak. While the relationship between NWS circulation and SSH variability warrants further investigation, it is beyond the scope of this paper.

The two spectral peaks (5–8 and 8–20 years) associated with the atmospheric variability (CtrlAtV) are similar to the peaks in the GC3.0 modelled NAO spectra (not shown), suggesting the NAO plays an important role in driving this variability. The NAO is a well-established driver to NWS SSH variability (Wakelin et al. 2003; Tsimplis et al. 2006; e.g., Su et al. 2014) and affects both the loading and the redistribution terms. The broad spectral peak (8–16 years) in SSH variability associated with the ocean boundary variance (CtrlOcV) is similar to the spectra of the GC3.0 modelled sub-polar gyre index (Hermanson et al. 2014) and modelled Atlantic Meridional Overturning Circulation (not shown). The sub polar gyre regulates the inflow of warmer and saline subtropical water into the North East Atlantic, leading to steric anomalies adjacent to the NWS (Chafik et al. 2019). The similarities of these spectra offer an intriguing insight into important drivers of NWS SSH variability. One of the benefits of the downscaling with the higher resolution, tidal, shelf seas model, is the improved representation of the shelf exchange processes, which may help modulate and control the propagation of the oceanic signals onto the shelf. Chafik et al. (2019) suggested that the along-shelf-break winds correlated with these exchange processes, and correlated with the NWS sea level variations. Both the remote drivers and the mechanisms of propagation of sea-level variations onto the shelf are important topics for future studies.

The final aim of this paper is to use the modelled variability to provide some context for coastal planners and managers using sea-level projections and estimates of sea-level variability to make decisions developing sea defence and adaptation schemes. We address this aim in two ways: quantifying the likely range of unforced sea-level change for a given period; and identifying the stretch of coast that a given tide-gauge based sea-level variability estimate is valid for.

Sea-level rise will not be as smooth as depicted in sea-level projections but will include variability on a range of timescales (e.g., Hinkel et al. 2019). This must be taken into consideration when using projections, as threshold levels may be exceeded sooner, and sea levels may be higher than projected, even if only for a period. We have quantified the distribution of trends, for a range of timescales, from our modelled annual mean SSH data. We choose the time periods from the modelled spectral peaks, rather than convenient round numbers. The largest changes in sea level occur on the 4-year timescale, when a shelf mean changes in sea level of 58.1 mm is possible (at the 95th percentile)—this is similar to the greatest likely global mean sea-level rise (50th percentile sea-level trend in RCP8.5 in the year 2100, ~ 12 mm/year, IPCC 2013). We find similar absolute changes for each time period considered, and so the rates vary considerably. If, for example, a 100-year time series included a 100 mm step change, the 5-, 10-, 20- or 50-year trend encompassing that step would have very different rates of changes, while having the same absolute change. This would also occur if the trends are driven by the occasional relatively large anomaly, especially when considering the tails of the distribution (i.e., the 95th percentile). The 21-year trend distribution is of a similar length to the altimeter record length. Sterlini et al. (2017) showed the altimeter observed North Sea linear sea level trend ranged from 1.3–3.9 mm/year from 1993–2014, with greatest levels off the German/Danish coasts (3.6 mm/year) and around Scotland (2.6 mm/year). Our results suggest that unforced variability may allow up to 2.0 mm/year and 1.8 mm/year in these regions respectively (Additional Material Fig. 3).

Traditionally, the best estimates of interannual to interdecadal variability come from long tide-gauge records. However, there are few old tide gauges, and so estimates are spatially sparse, with large lengths of the coast where no estimate can be made. There is a relatively high level of coherence around the coast at the interannual timescale, with a typical coastal grid box being representative of the adjacent coast. SSH variability of most coastal grid boxes can be used to represents 215 km of the coast of Great Britain, and 399 km of the European Coast (the respective median values for each coast). The more complex coastline of Great Britain tends to reduce the length scales, although in locations such as Newlyn (505 km) this complexity allows a longer stretch of coast to occur within smaller area. Brest and Santander represent the longest stretches of coast (887 km and 841 km respectively), both of which benefit from being on the relatively open coasts, while Brest (like Newlyn) is aided by being on the end of a peninsula. Further northeast along the European coast, the tide gauges of the southern North Sea tend to represent a smaller length of coast, which reflects how the amplitude of variability varies spatially in this region. We have focused on the redistribution term, as when using the total sea level (loading + redistribution) the spatially homogeneous loading component increases these coastal lengths considerably—we do include these values in the Table 12 and Table 13 in additional materials. We also note that filtering out the higher frequencies does not really change these coastal lengths while using the redistribution term but increase them when using the total sea level. This analysis could be extended to consider the length and quality of the tide gauge record—currently an estimate from a 100-year tide gauge record has the same “weighting” as a 15-year tide gauge. The current accepted approach is to estimate the local sea-level variance from a local long tide-gauge record. Our results support this methodology and suggest the portion of coast that a given tide gauge may be used.

Several methodological simplification and assumptions have led to caveats in this study. We do not consider variability in the Baltic exchange boundary condition, or in the riverine forcings, both of which are climatological in this study. We do not run additional barotropic model simulations which would allow us to separate the remote steric term from the local bottom pressure term (the method used by Chen et al. 2014). Both our models are Boussinesq, and so the steric effects on sea level must be diagnosed from rather than prognosed by the model. We do not consider the global mean thermosteric sea-level variability or trend—this is an important point for sea-level projections (Hermans et al. 2020b), but is less important for the sea-level variability in the climate control simulations used here. We focus on the interannual timescale as the baroclinic shelf seas model has not been optimised to simulate surge events – if we were interested in the higher frequency events, a tuned barotropic surge model may be more appropriate (e.g., O’Neill et al. 2016).

6 Conclusions

We conclude that:

• Our modelling system is able to simulate the inter-annual to decadal sea level variability on the NWS.

• Most of the NWS interannual SSH variability is forced by the interannual variability of the atmospheric forcings (CtrlAtV).

• SSH variability associated with ocean boundary variance (CtrlOcV) is coherent and so has no spatial pattern on the NWS. It is roughly half the magnitude of SSH variance associated with the atmospheric forcings (when averaged over the NWS).

• There is a spatial pattern to the SSH variability, with greatest variance in the eastern North Sea and around Scotland. This is associated with the variance in the atmospheric forcings.

• NWS SSH variability can be broken down into a shelf loading and redistribution term, and the covariance of the two. When averaging over the shelf, the redistribution term is roughly a third that of loading term for the atmosphere driven SSH variance (CtrlAtV) and is near zero for the ocean driven SSH variance (CtrlOcV). Therefore, when averaged over the shelf, the Ctrl redistribution term is roughly 20% of loading term.

• The covariance term is an important component of the spatial pattern, and the variance around Scotland is predominantly associated with covariance, perhaps reflecting its co-location with the main pathway of water into the North Sea.

• These terms are frequency dependent. The atmosphere boundary driven variance tends to be greater at high frequency. The correlations between the shelf loading and redistribution terms appear to change with frequency, and therefore, so does the balance between the redistribution term, and covariance term, as shown in the eastern North Sea.

• Unforced variability can lead to short terms changes in sea level, depending on location and timescale. The largest likely (95th percentile) modelled change occurs on the 4-year timescale, with a possible shelf mean increase of ~ 60 mm over 4 years. This is greatest in the eastern North Sea, where up to 90 mm is modelled. When considering sea-level projections, unforced sea-level variability, as quantified here, must be considered.

• The simulated tide gauges are typically coherent along an adjacent stretch of coast. The median length of this coast (with correlations greater than 95% and standard deviations within 15%) is 215 km for Great Britain, and 399 km for the (modelled) European coastline. Including the (spatially uniform) loading term increases these distances substantially. In the real world, other processes not modelled, may reduce this coherence.

Appendix

1.1 The non-linear free surface and the inverse barometer effect

Dynamic sea level (DSL) is the regional departure from the global-mean sea level, associated with changes in the ocean density and circulation (Roberts et al. 2016; Gregory et al. 2019), and is the variable stored as “zos” in CMIP5 and CMIP6. In this study, the surface elevation simulated by our models is equivalent to DSL, which we refer to as Sea Surface Height (SSH). By focusing on DSL, we do not consider the contributions from variability in the global mean thermosteric sea level, vertical land motion, change in the earth’s geoid and mass exchange with the cryosphere (Roberts et al. 2016).

Both the NEMO ocean component of GC3.0 and CO6 use a non-linear free surface. In the free surface formulation within NEMO, the variable η_B describes the shape of the air-sea interface (using the Boussinesq approximation). This prognostic variable is the vertical average of the kinematic surface condition:

$$\frac{\partial {\eta }_{B}}{\partial t}=-\nabla \bullet \left[\left(H+{\eta }_{B}\right)\stackrel{-}{{U}_{h}}\right] +\frac{{Q}_{m}}{{\rho }_{0}}$$

(6)

where $\nabla$ is the divergence operator (H is depth, U_h is the depth mean current), Q_m is the mass flux and ρ₀ is the reference density (set to 1020 kg m⁻³ in CO6 to reflect the typical shelf densities). This equation allows barotropic external gravity waves (such as tides). These waves have quite high phase speeds, so their timescales are short relative to most other modelled processes. In order to represent them efficiently, time-splitting is used. As the sea surface height changes, the entire vertical grid expands and contracts (Madec and the NEMO Team 2016).

Both models use the Boussinesq approximation (see Sect. 2.1.3) which means that steric expansion is not modelled prognostically but must be diagnosed (see Sect. 11.2). However, local changes in density (the local steric effect, see Sect. 2.1.3) impact sea level through changes in pressure.

In CO6, the inverse barometer effect (Stammer and Hüttemann 2008) η_ib, is included within the model calculation, and is calculated as

$${\eta }_{ib}=-\left(\frac{{p}_{s}-101000}{g{\rho }_{0}}\right)$$

(7)

where p_s is the surface pressure field (with a reference value of 101,000 Nm⁻²), g is gravity and ρ₀ is the reference density. η_ib is subtracted from internally modelled η and from the η external boundary conditions. The inverse barometer effect is not modelled in GC3.0.

1.2 Diagnosis of local steric expansion and bottom pressure changes

The local steric height and bottom pressure are diagnosed within CO6 at every model time step. Local steric height is (η_s) is calculated as:

$${\eta }_{s}=-\underset{H}{\overset{\eta }{\int }}\left(\frac{\rho \left(T,S,p\right)-{\rho }_{0}}{{\rho }_{0}}\right)dz=-\frac{1}{{\rho }_{0}}\left(\underset{H}{\overset{\eta }{\int }}\rho \left(T,S,p\right)dz- \underset{H}{\overset{\eta }{\int }}{\rho }_{0}dz\right)$$

(8)

where ρ and ρ₀ are the in-situ density and reference density, η the local sea surface height and H the local depth. It effectively shows the difference in water column height (H + η) between the in-situ density and reference density if the mass of the water column remained constant.

The bottom pressure (p_z) is calculated as:

$${p}_{z}= { g\underset{H}{\overset{\eta }{\int }}\rho \left(T,S,p\right)dz= \rho }_{0}g\left(\eta +H+{\eta }_{s}\right){10}^{-4}$$

(9)

where g is gravity. This effectively integrates the in-situ density over the water column.

1.3 Fraction of variance

We estimate the fraction-of-variance (FOV) to explain the relative importance of local bottom pressure and local steric changes (defined in Sect. 11.2) in the sea surface height. Roberts et al. (2016), their Eq. (5) defined the fraction of variance as

$${F}_{process}=1-\frac{{{var}} \, \left(ssh-{ssh}_{process}\right)}{{var} \, \left(ssh\right)}$$

(10)

where ssh is the SSH time series at each grid box. ssh_process is the component of SSH associated with either changes in local bottom pressure or local steric changes. F_process is the FOV associated with local bottom pressure and local steric changes respectively.

Recalling the rules of addition and subtraction of variances:

$${{var}}\left(X+Y\right)= {{var}}\left(X\right)+{{var}}\left(Y\right)+2{{cov}}\left(X,Y\right)$$

(11)

$${{var}}\left(X-Y\right)= {{var}}\left(X\right)+{{var}}\left(Y\right)-2{{cov}}\left(X,Y\right)$$

(12)

given (12), (10) may be rewritten as

$$\begin{aligned} F_{process} & = 1 - \frac{{{{var}}\left( {ssh} \right) + {{var}}\left( {ssh_{process} } \right) - 2{{cov}}\left( {ssh,ssh_{process} } \right)}}{{{{var}}\left( {ssh} \right)}} \\ & = 1 - \left( {\frac{{{{var}}\left( {ssh} \right)}}{{{{var}}\left( {ssh} \right)}} + \frac{{{{var}}\left( {ssh_{process} } \right)}}{{{{var}}\left( {ssh} \right)}} - \frac{{2{{cov}}\left( {ssh,ssh_{process} } \right)}}{{{{var}}\left( {ssh} \right)}}} \right) \\ & = \frac{{2{{cov}}\left( {ssh,ssh_{process} } \right)}}{{{var}\left( {ssh} \right)}} - \frac{{{{var}}\left( {ssh_{process} } \right)}}{{{{var}}\left( {ssh} \right)}} \\ \end{aligned}$$

(13)

F > 0 when 2cov(ssh,ssh_process) > var(ssh_process). When the 2cov(ssh,ssh_process)—var(ssh_process) = var(ssh), F = 1.

Here, we consider ssh to be the sum of two processes, so the sum of the two fractions of variance can be expressed as

$${F}_{1}+{F}_{2}=\frac{2( {{cov}}\left(ssh,{ssh}_{1}) +{{cov}}(ssh, {ssh}_{2}\right))}{{{var}}\left(ssh\right)}-\left(\frac{{{var}}\left({ssh}_{1}\right)+{{var}}\left({ssh}_{2}\right)}{{{var}}\left(ssh\right)}\right)$$

(14)

Given the additive law of covariances

$${{cov}}\left(A,B\right)+{{cov}}(A,C)= {{cov}}\left(A,B+C\right)$$

(15)

Equation (14) may be written as

$${F}_{1}+{F}_{2}=\frac{2( {{cov}}\left(ssh,{ssh}_{1}+{ssh}_{2}\right))}{{{var}}\left(ssh\right)}-\left(\frac{{{var}}\left({ssh}_{1}\right)+{{var}}\left({ssh}_{2}\right)}{{{var}}\left(ssh\right)}\right)$$

(16)

When ssh = ssh₁ + ssh₂ and ssh₁ and ssh₂ are uncorrelated, for example represented by Guassian distrubuted random number N(μ, σ), with μ mean and σ standard diviation:

$$\begin{gathered} ssh = ssh_{1} + ssh_{2} \hfill \\ ssh_{1} = {\text{N}}\left( {0,{1}} \right) \hfill \\ ssh_{2} = {\text{N}}\left( {0,{1}} \right) \hfill \\ \end{gathered}$$

(17)

cov(ssh₁, ssh₂) = 0 (as they are uncorrelated). Given (15) and (17), cov(ssh, ssh₁ +  ssh₂) = cov(ssh, ssh) = var(ssh). Given (17) and (11), var(ssh₁) + var(ssh₂) = var(ssh), and so (16) reduces to 2 − 1 = 1. Therefore, for two uncorrelated ssh processes that add exactly to give the total sea-level, the two fraction of variances adds to 1.

However, if the processes are correlated, or there is a third process (or noise), the sum of the fraction may not equal zero. If

$$\begin{gathered} ssh = ssh_{1} + ssh_{2} \hfill \\ ssh_{1} = {-}ssh_{2} = {\text{N}}\left( {0,{1}} \right), \hfill \\ \end{gathered}$$

(18)

the variance of ssh = 0, so (16) tends to -infinty. When you add noise to ssh (with an amplitude 10% of ssh₁), i.e.,

$$\begin{gathered} ssh_{1} = {-}ssh_{2} = {\text{ N}}\left( {0,{1}} \right) \hfill \\ ssh = ssh_{1} + ssh_{2} + {\text{N}}\left( {0,0.{1}} \right) = {\text{N}}\left( {0,0.{1}} \right) \hfill \\ \end{gathered}$$

(19)

the var(ssh) is √10% of var(ssh₁) and var(ssh₂), so the second term of (16) become 200 ((1 + 1)/√0.1). As ssh is now uncorrelated with ssh₁ +  ssh₂ (as it is only the guassian noise), the cov(ssh, ssh₁ +  ssh₂) is 0, so the sum of the fraction of variances in (16) the reduces to 0 − 200 = − 200.