Twenty-first century marine climate projections for the NW European shelf seas based on a perturbed parameter ensemble

The northwest European shelf (NWS) seas are environmentally and economically important, and an understanding of how their climate may change helps with their management. However, as the NWS seas are poorly represented in global climate models, a common approach is to dynamically downscale with an appropriate shelf sea model. We develop a set of physical marine climate projections for the NWS. We dynamically downscale 12 members of the HadGEM3-GC3.05 perturbed parameter ensemble (approximately 70 km horizontal resolution over Europe), developed for UKCP18, using the shelf sea model NEMO CO9 (7 km horizontal resolution). These are run under the RCP8.5 high-greenhouse-gas-emission scenario as continuous simulations over the period 1990–2098. We evaluate the simulations against observations in terms of tides, sea surface temperature (SST), surface and near-bed temperature and salinity, and sea surface height. These simulations represent the state of the art for NWS marine projections. We project an SST rise of 3.11 °C (± 2σ = 0.98 °C) and a sea surface salinity (SSS) freshening of −1.01 (± 2σ = 0.93; on the (unitless) practical salinity scale) for 2079–2098 relative to 2000–2019, averaged over the NWS (approximately bounded by the 200 m isobar and excluding the Norwegian Trench, the Skagerrak and Kattegat), a substantial seasonal stratification increase (23 d over the NWS seas), and a general weakening of the NWS residual circulation. While the patterns of NWS changes are similar to our previous projections, there is a greater warming and freshening that could reflect the change from the A1B emissions scenario to the RCP8.5 concentrations pathway or the higher climate sensitivity exhibited by HadGEM3-GC3.05. Off the shelf, south of Iceland, there is limited warming, consistent with a reduction in the Atlantic Meridional Overturning Circulation and associated northward heat transport. These projections have been publicly released, along with a consistent 200-year present-day control simulation, to provide an evidence base for climate change assessments and to facilitate climate impact studies. For example, we illustrate how the two products can be used to estimate climate trends, unforced variability and the time of emergence (ToE) of the climate signals. We calculate the average NWS SST ToE to be 2034 (with an 8-year range) and 2046 (with a 33-year range) for SSS. We also discuss how these projections can be used to describe NWS conditions under 2 and 4 °C global mean warming (compared with 1850–1900), as a policy-relevant exemplar use case.

1 Introduction

The northwest European shelf (NWS) seas are economically and environmentally important (Pugh, 2008; Tinker et al., 2018). They are quasi-isolated from, and behave differently to, the North Atlantic, and they have a different response to climate change (Holt et al., 2010; Wakelin et al., 2009). They are generally poorly represented in global climate models (GCMs) and global ocean models due to limited vertical and horizontal resolution and the absence of key processes, such as tides (Tinker et al., 2022). Therefore, assessing the NWS response to climate change directly from GCMs may be inappropriate for many applications (Hermans et al., 2020b; Tinker et al., 2022; Mathis et al., 2013). A common approach to mitigate these shortcomings is to dynamically downscale GCM climate projections with a higher-resolution shelf sea model, which includes the additional relevant processes (e.g. Holt et al., 2010; Tinker et al., 2016, 2020).

The marine component (Palmer et al., 2018) of the recent national UK Climate Projections, UKCP18 (Murphy et al., 2018), had a sea level focus (see their table 5.1). UKCP18 provided mean sea level projections (for the 21st century and exploratory extended projections to 2300; Palmer et al., 2020), estimates of the change in the surge and wave climate (Howard et al., 2019), and a quantification of the present-day sea level variability (Tinker et al., 2020). The climate projections of ocean water properties for the NWS from UKCP09 (Lowe et al., 2009) and the Minerva projections (Tinker et al., 2015, 2016) were not updated in UKCP18. The UK's third Climate Change Risk Assessment (CCRA3) used UKCP18 as a primary source of evidence. However, as the NWS climate projections had not been updated, assessment of those aspects was based on the older UKCP09/Minerva NWS projections. As well as being based on outdated science, use of the earlier NWS projections may introduce inconsistencies with the wider CCRA3 findings based on UKCP18. With the evidence call for the upcoming fourth CCRA approaching, we have developed a set of NWS marine climate projections, based on and consistent with the UKCP18 marine projections.

Beyond the UK, there are numerous legislative mechanisms that support the protection and management of the NWS marine environment. At the European level, examples include the following: the EU Common Fisheries Policy (CFP, 2013); the EU Habitats Directive (1992); the Oslo–Paris (OSPAR) “Convention for the Protection of the Marine Environment of the North-East Atlantic” (1992; see OSPAR, 2010); and the EU Marine Strategic Framework Directive (MSFD, 2008). There are also important international treaties and conventions such as the United Nations Framework Convention on Climate Change (UNFCCC, 1992) and the Convention on Biological Diversity (CBD, 1992). These are all implemented at the UK level (Frost et al., 2016), and all have varying and overlapping goals and targets linked to the protection and monitoring of the marine environment. While not many of these policies explicitly mention climate change, they capture its effect through natural variability and broader environmental change (Frost et al., 2016). For example, the MSFD aims to achieve “Good Environmental Status” (GES), which varies by region across the EU waters. After an initial series of surveys to define GES, the MSFD has driven a series of monitoring programmes to chart the progress towards GES and the implementation of a set of management strategies (programmes of measures) to help its achievement. Marine climate projections help inform all three steps in this process; without an idea of how the marine climate may change, it is difficult to make informed management decisions, and monitoring must consider the background climate change to help interpret any observations. There is also an increased focus on developing conservation and planning legislation that takes into account future projections of climate change (Queirós et al., 2021; MCCIP, 2023). Effective marine management and adaptation planning requires information on the past and present environment through monitoring and information on future changes from climate projections.

In this paper (and dataset), we downscale the UKCP18 twelve-member perturbed parameter ensemble (PPE), based on HadGEM3-GC3.05, run under the RCP8.5 high-emissions climate change scenario. The PPE simulations are downscaled with the shelf seas version of NEMO 4.04 (Coastal Ocean version 9, CO9) as transient simulations for the period 1990–2098. We use the methodology and evaluation developed by Tinker et al. (2020) to produce a set of projections consistent with their estimate of NWS unforced variability. These updated NWS climate projections (hereinafter NWSPPE) are released, with the present-day control simulation (hereinafter PDCtrl) of Tinker et al. (2020), in time for use in the UK's fourth CCRA.

1.1 Choice of RCP8.5

Our set of projections are based on the RCP8.5 concentration trajectory. RCP8.5 is a relatively high-impact business-as-usual scenario from the CMIP5 suite of models rather than the more recent CMIP6, which is based on the Shared Socioeconomic Pathways (SSPs). RCP8.5 has very similar total radiative forcings to the SSP5-8.5 scenario (Tebaldi et al., 2021), although RCP8.5 has a slightly weaker global temperature response, attributed to a lower CO₂ concentration (Fyfe et al., 2021). The choice of forcing scenario is motivated primarily by the desire for a set of NWS marine climate projections consistent with the latest set of UK Climate Projections (UKCP18), which were run under RCP8.5 and before the SSP scenarios became available. This consistency allows researchers to look across both the land and marine domains in multi-variate space in a way that has previously not been possible to facilitate, e.g. consideration of compound hazards and the combined effects of multiple climate-impact drivers.

Recent studies have criticised the use of RCP8.5 in climate projections, particularly in terms of its apparent low likelihood, given emission reduction pledges associated with the Paris Agreement (Hausfather and Peters, 2020). However, there are several scientific reasons why this remains a useful scenario in the context of policy-relevant climate information. Firstly, it has a high signal-to-noise ratio and can therefore better separate the forced climate response from internal variability. Secondly, RCP8.5 can be readily translated into warming levels, which was the primary basis of the last UK Climate Change Risk Assessment, and appears prominently in the IPCC Sixth Assessment Report (AR6; IPCC, 2021a). Thirdly, risk-based decision-making requires a comprehensive picture of the future risk landscape, including higher warming levels associated with any combination of emissions back-tracking, positive carbon-cycle feedbacks and high climate sensitivity (IPCC, 2021b). Finally, high-emission scenarios such as RCP8.5 provide a useful baseline scenario from which the benefits of mitigation action and avoided costs can be assessed.

2 Models and methods

We use CO9 to downscale the UKCP18 HadGEM3-GC3.05 PPE (Yamazaki et al., 2021). This ensemble is downscaled as transient simulations from January 1980–November 2099, but the first 10 years are considered spin-up and so are discarded. This leads to a 12-member PPE downscaled for the NWS (NWSPPE), from January 1990 to December 2098, for RCP8.5.

The NWSPPE is consistent with PDCtrl, which can provide an estimate of a different source of uncertainty, and so the two datasets are released together.

Here we give a brief overview of HadGEM3-GC3.05 (our driving GCM) and NEMO 4.0.4 CO9 (our regional ocean model), its treatment of tides, and where the model code is available. A detailed description of HadGEM3-GC3.05 and its evaluation for the NWS seas is given in Tinker et al. (2020). We then describe how HadGEM3-GC3.05 has been run as the perturbed parameter ensemble before describing some of our methodologies.

2.1 HadGEM3-GC3.05

HadGEM3-GC3.05 (Williams et al., 2018; Yamazaki et al., 2021) is the Met Office GCM used in the 2018 UK Climate Projections (UKCP18) and very similar to the model submitted to CMIP6 (Coupled Model Intercomparison Project Phase 6; Eyring et al., 2016). HadGEM3-GC3.05 is based on the Met Office Unified Model atmosphere (Walters et al., 2019) (with approximately 70 km resolution of Europe) and the ocean model NEMO, run on the ORCA025 grid (∼ $\mathrm{1}/\mathrm{4}$° horizontal resolution on a tri-polar grid). HadGEM3-GC3.05 uses the CICE sea ice model (Hunke et al., 2015) with the GSI8.0 sea ice configuration (Ridley et al., 2018) on the same ORCA025 grid. HadGEM3-GC3.05 does not simulate tides or other important shelf sea processes. HadGEM3-GC3.05 is run with a 360 d calendar, with 12 months of 30 d. This has little impact on the climate projections but is not directly compatible with observed tides, which are constrained by astronomical frequencies. This is discussed later.

We use HadGEM3-GC3.05 model output from the atmosphere and the ocean model to drive our downscaled shelf sea simulations. HadGEM3-GC3.05 atmospheric surface fluxes of heat, fresh water and momentum are used, as are ocean temperature, salinity, currents, and sea surface height from the HadGEM3-GC3.05 ocean model. Tinker et al. (2020) lists all the HadGEM3-GC3.05 variables used to drive our simulations, and their frequency, in their Table S9.

The HadGEM3-GC3.05 perturbed parameter ensemble

HadGEM3-GC3.05 was run as a perturbed parameter ensemble (PPE) for the UKCP18 climate projections. To explore the uncertainty associated with the choice of parameters and parameterisations within the GCM, a PPE perturbs the parameters (within a reasonable range) and chooses different parameterisations for each ensemble member. As more ensemble members are added, the sensitivity to these parameters increases the ensemble spread. This spread was designed to span the climate uncertainty associated with these uncertain parameters. The GCM simulations used and the PPE framework are described in Sexton et al. (2021) and Yamazaki et al. (2021).

Yamazaki et al. (2021) ran a 25-member HadGEM3-GC3.05-PPE for the period 1900–2100, using CMIP5 historical and RCP8.5 emissions. The variants had different combinations of perturbations to 47 parameters in the model's atmosphere, land and aerosol schemes, and they were selected by running a set of cheap, coarser-resolution atmosphere-only simulations from a large sample of nearly 3000 variants. Poor-performing variants were filtered out by assessing retrospective 5 d weather forecasts and simulations of 2004–2009 with prescribed sea surface temperature (SST) and sea ice. Further idealised atmosphere-only experiments were then used to select 25 variants which were as diverse as possible in terms of their atmospheric climate feedbacks and their regional responses to warmer SSTs (Sexton et al., 2021). These variants were then run with HadGEM-GC3.05 for historical and future simulations from 1900–2100 under RCP8.5. Thus, 10 members were dropped as being too cold by 1970 or when evaluated against 13 CMIP5 models over the historical period (1900–2005). The 15 resulting ensemble members gave plausible yet diverse atmosphere and ocean model behaviours (Yamazaki et al., 2021), and they were used within the UKCP18 climate projections. Of this 15-member ensemble, 3 ensemble members exhibited unrealistic Atlantic Meridional Overturning Circulation (AMOC) behaviour and so were downscaled for the NWS but were excluded from our analysis. The 12 remaining downscaled PPEs are referred to as the NWSPPE.

Heat and salinity flux adjustments were applied to the HadGEM3-GC3.05 PPE to prevent the surface climate from drifting too far from the realistic state, despite the parameter perturbations giving rise to a range of net top-of-atmosphere radiative fluxes. Details of how the flux adjustments were applied are described in Yamazaki et al. (2021).

A recent study compared UKCP18 PPE simulations for RCP8.5 to a selection of available CMIP6 models (Murphy et al., 2023). The evaluation of the climatology, large-scale modes of variability and the historical change in global mean surface temperature showed that the overall level of performance from the UKCP18 global PPE simulations was competitive with the CMIP6 ensemble. Both sets of simulations exhibit a range of biases across their individual members, with evidence of a systematic component to the biases for some variables (Murphy et al., 2023). When looking at the multidecadal average changes, the CMIP6 ensemble develops a broad range of warming signals in response to SSP5-8.5 (e.g. 2061–2080 relative to 1981–2000), commensurate with a broad range of climate sensitivity values among models (e.g. IPCC, 2021c). This wide CMIP6 range shows substantial overlap with those of the UKCP18 global PPE for both surface temperature and precipitation, confirming the continued importance of accounting for the relevant uncertainties in impact studies. Both UKCP18 global PPE and CMIP6 estimates span positive and negative changes in winter and summer North Atlantic Oscillation (NAO), although the UKCP18 shows a slight narrowing of the range of summer NAO compared with CMIP6. The study highlights the importance of assessing sources of uncertainty (including structural modelling uncertainty), which is often overlooked in regional marine climate projections (see Sect. 6 for further discussion).

2.2 NEMO AMM7

We use the shelf seas version of NEMO (O'Dea et al., 2017) to dynamically downscale HadGEM3-GC3.05. It is run on the AMM7 NWS domain, which extends from 40°4^′ N, 19° W to 65° N, 13° E with a 7 km horizontal resolution with hybrid terrain-following vertical levels (s-levels; Siddorn and Furner (2013) for PDCtrl and Bruciaferri et al. (2018) and Wise et al. (2022) for the transient projections).

Two slightly different versions of the model are run for the present-day control simulation and the climate projections. Coastal Ocean version 6 (CO6), based on NEMO version 3.6, was used in the PDCtrl – see Tinker et al. (2020) for full details. Coastal Ocean version 9 (CO9, NEMO 4.04) is used for the NWSPPE climate projections. Tests showed little difference between the two models, so we consider the two versions to be consistent. We use “COx” when we refer to both model versions.

COx is driven by several different classes of boundary conditions. The atmospheric and ocean boundary conditions are provided by HadGEM3-GC3.05 and are interpolated onto the COx grid. The Baltic exchange and the rivers are specified by observation-based climatologies. Full details of how these are implemented and their frequencies are given in Tinker et al. (2020) and their Table S9 in particular.

2.2.1 NEMO AMM7 tides

The standard version of NEMO introduces errors into the tides when using the 360 d climate. It resets the tides at the beginning of each month and then runs sequentially. This leads to a jump at the end of the months with lengths other than 30 d. This does not seem to cause the model to crash but must affect the internal model solution. It also, on average, changes the M₂ period and so rendered standard tidal analysis useless.

We have adapted the tidal module within NEMO to work with the 360 d calendar, with three different methods: “reset”, the standard method, where the tides are reset to every month; “drift” (which is used here), where we allow the tides to run continuously from a given date but the annual tides (such as the equinoctial tides) drift by about 5 d every year compared to the Gregorian calendar; and “compress”, where we alter one of the astronomical frequencies within the tide module – this changes the M₂ (and most other) tidal frequency by ∼ 2 s but keeps the equinoctial tides with their correct timings. We describe these three methods and evaluate the impact they have on the tidal behaviour in Appendix B. In these simulations, we use the drift method, where the model time, in number of seconds since 1 January 1980 (in the 360 d calendar), is counted forward with the Gregorian calendar to give the tidal conditions.

Both model versions use the same tidal constituents and tidal boundary conditions; however, in the later version of NEMO (4.04) the default tidal Love number was changed from 0.7 to 0, effectively turning off the internal tide-generating potential. This was not realised until the projections had been completed. In a small domain like the NWS, most tidal energy is propagated in from the boundary, and very little energy is internally generated; however, this slightly affects the tidal range. A sensitivity test was run (running the unperturbed ensemble member with the tide-generating potential on and off), showing that this only had a small impact on the tidal range (Fig. S1), changing the amplitude of the M₂ tide by up to 4 cm (which is up to about 3.5 m). This was not considered significant for the climate projections, although in a different context (such as operational surge forecast systems) this could be more of a problem.

2.2.2 NEMO COx GitHub configurations

Both NEMO configurations used in this study are available on GitHub.

NEMO version 3.6 CO6 (AMM7)

https://github.com/hadjt/NEMO_3.6_CO6_shelf_climate (last access: 20 May 2024)

Development of this configuration has frozen. The latest version available on GitHub was used in PDCtrl.

NEMO version 4.0.4 CO9 (AMM7 and AMM15)

https://github.com/hadjt/NEMO_4.0.4_CO9_shelf_climate/ (last access: 20 May 2024)

Development of this configuration has continued since the NWSPPE was run.

The version used for NWSPPE is available:

https://github.com/hadjt/NEMO_4.0.4_CO9_shelf_climate/tree/bddd0f68980632229c5afef9772c9fd0d0d6e930 (last access: 20 May 2024)

A fix has since been added to set the default tidal Love number to 0.7:

https://github.com/hadjt/NEMO_4.0.4_CO9_shelf_climate/tree/4ef0ec5e0c20a9aa88b42f395cc9b1bfd689a221 (last access: 20 May 2024)

2.3 Present-day control simulation

HadGEM3-GC3.05 was also run as a present-day control simulation for UKCP18 and underpinned Tinker et al. (2020). It was run for 270 years while simulating a stable climate equivalent to the year 2000. This was forced with greenhouse gases and aerosols from the year 2000 (as fixed values or annually repeating seasonal cycles). This allowed the model to simulate the natural, unforced climate variability that would have occurred during the year 2000 and the absence of climate change. This simulation should capture the range of climate variability associated with different phases of climate modes such as the El Niño–Southern Oscillation, North Atlantic Oscillation and Atlantic Multidecadal Variability. The first 70 years of the simulation still show trends associated with the deep ocean spinning up, so the last 200 years were used for analysis by Tinker et al. (2020) and form the second part of this dataset. We refer to the downscaled control run as PDCtrl.

2.4 Copernicus Marine NWS reanalysis (RAN)

A Met Office-led consortium provides a marine reanalysis for the NWS to the Copernicus Marine Service (product number: NWSHELF_MULTIYEAR_PHY_004_009). This ran from 1993 to the present, using the same regional shelf sea model and domain as used here (NEMO CO6). By running a shelf sea model and assimilating observations (SST, temperature and salinity profiles and satellite altimeters), it can be considered a best-guess estimate of the current state of the NWS seas.

We make limited use of this NWS reanalysis (hereinafter RAN) but comparing its regional mean time series (Tinker et al., 2019) to those of our NWSPPE. This gives some early-century context to the NWSPPE temporal evolution in terms of interannual variability and trends.

2.5 Circulation, currents and transport

Residual currents are the small remaining currents after the large oscillatory tidal currents are averaged out. They are an important feature of the NWS seas and are an important determinant of the spatial distribution of heat, salt and nutrients. We use two methods to analyse the NWS residual circulation: spatial maps of barotropic (depth-mean) residual currents and time series of transport through pre-defined cross-sections. Following Tinker et al. (2022), we use 30 d mean barotropic currents (all simulated months have 30 d with a 360 d calendars) as a proxy for the residual circulation because they effectively average out the S₂ and M₂ tide – 720 h is 57.97 cycles of the M₂ tide (Tinker et al., 2022). We initially consider the spatial configuration of the NWS residual circulation using the magnitudes of the mean currents.

We also assess how much the currents vary about their mean using the residual current uncertainty ellipse method of Tinker et al. (2022). Normally distributed U and V currents together have a bivariate normal distribution (Fig. S2a), which follows a Gaussian surface rather than a Gaussian curve. Contours (of equal probability) on that surface enclose a given proportion of the data – the ellipse associated with 2.45 standard deviations encloses 95 % of the data (Fig. S2b and c). By fitting a bivariate normal distribution and these ellipses to the residual currents, we can use the ellipse properties to analyse the circulation in novel ways (Fig. S2c–f).

We can use the ellipse properties to compare the mean of the residual currents to their variability. When the mean is greater than their variability, the residual currents are relatively steady (not stopping or reversing); when the variability is greater than the mean, the residual current may have a preferential direction and strength, but at any given time, they could be going in any direction. We can also use this approach to compare two sets of residual current distributions from different times or from different ensemble members, and we see how similar they are by assessing the overlap (OVL) of their bivariate normal distribution (Fig. S2h and i). The OVL method finds the volume (or area) under two Gaussian surfaces (or curves) to quantify how similar they are. The volume under a single Gaussian function integrates to 1, so the integral under two Gaussian functions must integrate to ≤ 1 (see Fig. S2i) but also must be greater than 0, as Gaussian functions extend to infinity. OVL ≈ 1 when the two distributions are very similar (in terms of their means, variances and covariance) and OVL ≈ 0 when they are very different (see Fig. S2i). We extend the OVL method of Tinker et al. (2022) to quantify how similar the residual current distributions are across the ensemble (ens_OVL). By finding the volume under the bivariate normal distribution for all 12 ensemble members, we get a pointwise measure of the similarity of the residual circulation across the NWSPPE and can then show how this changes over the NWS, as well as over the 21st century.

Tinker et al. (2022) give a full description of the methodology and provide a Python toolbox on GitHub. We have included their figure explaining their methodology (Fig. S2).

Time series of transport through pre-defined cross-sections are calculated every time step at runtime. Thus, 76 cross-sections across the domain quantify the net, positive and negative transports of volume, heat and salt. We also output spatial maps of monthly mean transport (the instantaneous product of each current component with the depth of water) and are able to calculate the net volume cross-section transport – this allows us to add additional cross-sections post hoc.

We also use these volume transport cross-sections to help evaluate the model, using the same cross-sections and observations as used in UKCP09 (Lowe et al., 2009; Holt et al., 2010).

2.6 Regions and regional mean statistics

Time series are useful to show temporal evolution and how an ensemble behaves. We use a pre-defined region mask to calculate time series of the regional means (and other statistics). We use the region mask of Wakelin et al. (2012) as it divides the NWS into regions that make geographic and oceanographic sense. We also include a “shelf” region, by combining several NWS regions (excluding the Atlantic, Norwegian Trench, Skagerrak/Kattegat and Armorican shelf regions). When considering the sparse EN4 datasets, we combine these regions into larger validation regions (Tinker et al., 2019) to improve the data coverage and sample size. While there is relatively high data coverage within the North Sea, other NWS regions are still relatively data poor, especially given the need to average out the differing phases of variability, to compare the model and observation “climate” rather than “weather”. For completeness, we include the EN4 analysis on the Wakelin et al. (2012) mask in the Supplement (Fig. S3).

2.7 Normalised bias

During our evaluation, we calculate the normalised bias between the model and observations to show where the observations sit within the distribution of the ensemble. To do this, we subtract the observations from the ensemble mean and divide by the ensemble standard deviation:

$$\begin{array}{}\text{(1)}& \text{NormBias }={\displaystyle \frac{\text{Ens\hspace{0.17em}\hspace{0.17em}Mean}-\text{Obs}}{\text{Ens\hspace{0.17em}\hspace{0.17em}Std}}}\end{array}$$

This gives the number of (ensemble) standard deviations the observations are above or below the ensemble mean. We typically consider a value with an absolute normalised bias less than 2 to be within the NWSPPE.

2.8 Potential energy anomaly

We use the potential energy anomaly (PEA; Simpson and Bowers, 1981) to quantify of the strength of the stratification. PEA is equivalent to the amount of energy required to fully mix a stratified water column.

PEA is quantified internally by NEMO:

$$\begin{array}{}\text{(2)}& \mathrm{PEA}=-{\displaystyle \frac{g}{h}}\underset{z=-h}{\overset{\mathrm{0}}{\int }}z\left(\mathit{\rho }(T,S)-\mathit{\rho }(\stackrel{\mathrm{‾}}{T},\stackrel{\mathrm{‾}}{S})\right)\mathrm{d}z,\end{array}$$

where h is the depth of water (limited to the upper 400 m), g is gravity, z the vertical coordinate (positive upwards), and an over bar represents the depth average.

The thermal and haline components of PEA (PEAT and PEAS) can also be calculated by Eq. (2) by replacing the S with $\stackrel{\mathrm{‾}}{S}$ for PEAS and T with $\stackrel{\mathrm{‾}}{T}$ for PEAT.

We follow the Holt et al. (2010) definition of stratification, where the water column is stratified when PEA > 10 J m⁻³ and (mixed layer depth) MLD < 50 m.

2.9 Stratification initialisation and duration

Daily mean fields of PEA and MLD are used to assess the timing and duration of stratification.

For each year, and each grid box, we identify stratified days (PEA > 10 J m⁻³, MLD < 50 m) and days where stratification initialises (a stratified day following a mixed day) and breaks down (a mixed day following a stratified day). We note the grid boxes that are stratified throughout the year and those that are mixed throughout the year. We then cycle through these periods of stratification and count their duration. We identify the longest single period of stratification and note the days of the year that it initialises and breaks down, as well as its duration. This discards short periods of stratification that may precede the main period of stratification or may occur after the main breakdown. We also note the number of periods of stratification in the year and the total days of stratification. This method assumes that seasonally stratifying regions are mixed in the winter; it does not account for stratified periods that continue from one year to the next. We then calculate 20-year means of the stratification duration and days of the year of the initialisation and breakdown. These data are then considered across the ensemble. This method is consistent with that of Holt et al. (2010) and Jardine et al. (2023).

The analysis was also repeated with a stratification threshold of the DFT (the DiFference in Temperature between the surface and near-bed temperatures, i.e. DFT = SST − NBT) being greater than 0.5 °C and PEAT > 10 J m⁻³, both of which gave a similar result.

2.10 Estimation of the time of emergence (ToE)

We estimate the time of emergence (ToE) of the climate signal for SST and sea surface salinity (SSS), using the method of Lyu et al. (2014). We compare a modelled estimate of the present-day unforced variability (the noise N), to the smoothed climate trend of an ensemble (the signal S) and ask when the ratio of the signal to noise is greater than a given threshold (2 standard deviations), if it does not drop below a threshold for the rest of the record and is at least 20 years long. We calculate the ToE for each ensemble member, and we then report the median value and the 16th–84th percentile range.

Table 1 Overview of observation datasets (and their analysis) used in the NWSPPE evaluation.

We estimate the noise as the present-day unforced variability by taking the standard deviation of the linearly detrended PDCtrl annual means. As the PDCtrl has a near-zero trend, detrending the PDCtrl has little effect on the ToE. We calculate the signal of each ensemble member by smoothing with a fourth-order polynomial and converting to anomalies by removing the 2000–2019 baseline period. We then find where the S:N ratio exceeds the threshold and remains above the threshold for the rest of the simulation. We make one further criterion that the ToE does not occur in the final 20 years of the projection (and so the emergence time is at least 20 years long), as they may not represent the true emergence of climate signals.

2.11 Salinity units

Both NEMO AMM7 and HadGEM3-GC3.05 use the 1980 equation of state (EOS-80; UNESCO, 1983), which works with the practical salinity scale (which has no units). In this study, we report salinities on the practical salinity scale.

3 Observational datasets

In this section, we describe the observational datasets used in our NWEPPE evaluation (Table 1). Tinker et al. (2020) evaluated PDCtrl, and we follow a similar approach here.

3.1 OSTIA

We use the OSTIA (Operational Sea-Surface Temperature and Sea-Ice Analysis) product (Roberts-Jones et al., 2012) to evaluate the NWSPPE SST. OSTIA is a largely satellite-based SST dataset on a horizontal grid of 0.25° which we bi-linearly interpolated onto the AMM7 model grid. The OSTIA assimilates several bias-corrected satellite products to reduce the bias of the overall product.

Figure 1 The (2000–2019) ensemble mean cotidal chart for the four largest constituents (a) M₂, (b) S₂, (c) N₂ and (d) K₂, as well as (e) the shallow water component (M₄) and (f) the offset (mean sea level). Amplitude (m) is represented by the colour map, with black contours matching the colour-bar tick labels. The grey contours give the phase in 45° intervals.

Figure 2 SST bias of the NWSPPE ensemble mean for 2000–2019 (for the four seasons) compared with the OSTIA SST. The hatching shows where the OSTIA SSTs are more than 2 (ensemble) standard deviations from the ensemble mean.

We calculate an SST climatology (mean and standard deviation) for 2000–2019 from monthly, seasonal and annual mean OSTIA data. This is compared with the NWSPPE ensemble mean for the equivalent period (Fig. 2). We also assess where the OSTIA 20-year mean sits within the NWSPPE (when processed equivalently), using the normalised bias in Eq. (1), asking how many NWSPPE ensemble standard deviations OSTIA is form the NWSPPE ensemble mean.

3.2 EN4

We evaluate the surface and near-bed temperature and salinity of the NWSPPE with the EN4 quality-controlled temperature and salinity profile dataset (Good et al., 2013). We use the individual observed vertical profiles of version EN.4.2.2.g10.

As the EN4 data are relatively sparse on the NWS, it is a complex dataset to use in a climate context, and so we follow a similar methodology to Tinker et al. (2020). There are not many EN4 locations where observations are available for every year for a given month and grid box, and so it is seldom possible to separate the observed climatology from the observed interannual variability. Therefore, rather than comparing the observed climatology from the distribution of NWSPPE climatologies (the distribution of the 20-year means from the NWSPPE), we compare the EN4 observation to the NWSPPE distribution of both interannual variability and ensemble variability (12 ensemble members and 20 years). This is a like-for-like analysis. This approach gives a reasonable spatial coverage of the NWS, albeit with a very low number of repeat samples.

EN4 evaluation methodology

This comparison is designed to evaluate the broad characteristics of any model biases. A more detailed comparison is beyond the scope of the present study.

Figure 3 Regional mean distribution of normalised model minus EN4 profile observations for SST (red), NBT (blue) and SSS (green). Each EN4 profile–model pair is used to calculate a normalised bias (Eq. 1, model minus ensemble mean divided by ensemble standard deviation), as shown in Figs. S4–S6. These data points are separated into distributions by validation regions (f) and month; the different regions are plotted as separate panels, and the months are separated along the x axis (showing first letter of the month name). Each distribution is plotted as the mean and ± 2 standard deviation, with two numbers – the number of points in each distribution is the upper number and the distribution mean is the lower number. The validation region mask is given in panel (f), which denotes (a) the North Sea, (b) the Celtic Seas (including Irish Sea, Celtic Sea and English Channel), (c) outer shelf region, (d) the Norwegian Trench (including the Skagerak and Kattegat) and (e) the oceanic regions. This figure is repeated with the Wakelin region mask in Fig. S3, giving greater granularity but smaller sample sizes.

Creating the NWSPPE distribution fields.

• 1.

  We create a NWSPPE climatology.

  • a.

    For each of the 12 months, we average all the monthly means from all the ensembles between 2000 and 2019.

  • b.

    As the EN4 dataset is sparse, we also calculate the ensemble standard deviation across all 240 points for each of the 12 months (20 years and 12 ensemble members).

Pre-processing the EN4 data.

• 2.

  We discard the profiles with QC flags equal to 4 (POSITION_QC, PROFILE_POTM_QC, PROFILE_PSAL_QC).

• 3.

  For each month and for each year, we do the following:

  • a.

    We assign each profile to the nearest COx grid box and linearly interpolate it onto the vertical s-levels.

  • b.

    We average the profiles if there are more than one for a given month.

  • c.

    We calculate the normalised bias as follows.

    $\mathrm{NormBias}=\frac{\left({\text{PPE\_mean}}_{i,j,z,m}-{\text{Obs}}_{i,j,z,y,m}\right)}{{\text{PPE\_std}}_{i,j,z,m}}$,

    where $i,j,z$ are the grid indices in three dimensions, and y and m are the year and month.

  • d.

    For each month and for each grid box, we add all the NormBias values from all the years, add the square of the NormBias and count the number of years with an observation (Sum_NormBias, SSq_NormBias, cnt_NormBias).

  • e.

    For each month, we calculate the mean and the standard deviation of the NormBias for each grid box (Mean_NormBias, Std_NormBias).

  • f.

    We combine Sum_NormBias, SSq_NormBias and cnt_NormBias for each month into seasons (DJF, MAM, JJA and SON) and annual means.

This gives us a gridded assessment (for each month, season and year) of where the EN4 observations sit within the NWSPPE (SST: Fig. S4; SSS: Fig. S5; NBT: Fig. S6). For each grid box with an EN4 observation, we can then give the number of NWSPPE standard deviations that it is from the ensemble mean. If the grid box (for a given season or month) has observations from different years, they are averaged. These sparsely sampled maps are hard to use quantitatively, and so regional summaries are also produced (Figs. 3 and S3).

3.3 ICES climate time series

While the EN4 data give a reasonable coverage of the NWS, very few points have enough observations to allow the climate to be quantified in terms of its mean and variability. We supplement the EN4 dataset with several long-time-series observations that expand our comparison of the NWSPPE simulations and observations.

The International Council for the Exploration of the Sea (ICES) observations support sustained-observation time series to allow for long-term monitoring of the ocean climate of the North Atlantic. These are formed by regular measurements of ocean temperature and salinity over decades (González-Pola et al., 2022). Furthermore, many ICES records are for the subsurface, which is not captured by satellites.

Table 2 ICES long-term climate time series used in the evaluation of the NWSPPE, giving the location, depth of water, and how the time series was calculated from the temperature and salinity depth profiles. (NAW: North Atlantic Water; MNAW: Modified North Atlantic Water; WCO: Western Channel Observatory; AFBI: Agri-Food and Biosciences Institute)

Figure 4 Comparison of the NWSPPE with long annual mean time series of in situ T and S from around the NWS (see text). The left column shows temperature (a–k), and the right column shows salinity. Only temperature was measured in panels (i–k). These records are a mix of surface, near-bed, water column mean (a–k) and the mean of the upper water column.

We evaluate the NWSPPE temperature and salinity against 11 ICES ocean climate time series (González-Pola et al., 2022; Hughes et al., 2018) listed in Table 2 (data retrieved from https://ocean.ices.dk/core/iroc, last access: 12 December 2023). Of these 11, 8 are from the NWS; 2 are from the north of the NWS, towards the Faroe Islands (Faroe Shetland NAW, Faroe Shetland MNAW), and 1 is from the Norwegian Trench (Utsira B). The time series are a mix of surface, near-bed, depth-averaged, and the average of a portion of the water column (or of a particular water body), and they are listed in Table 2. We therefore extracted the annual mean temperature and salinity water column data from the nearest model grid box, and we applied the equivalent depth processing. We compare data from 2000–2019 and consider whether the observations are within the range of the NWSPPE, as well as comparing the mean and the interannual variability (Fig. 4, Table S1). We calculate the normalised bias following Eq. (1) to show how much warmer/saltier the model is from the observations, in ensemble standard deviations. To assess the variability of the model and observations, we calculate the standard deviation of the ICES time series and for each of the 12 ensemble members, and then we calculate the normalised bias following Eq. (1).

3.4 Satellite sea surface height (SSH) data

We evaluate the NWSPPE SSH following the methodology of Tinker et al. (2020), where we refer the reader for evaluation of PDCtrl SSH. We compare the mean pattern of the modelled sea level to satellite altimetry mean dynamic topography (MDT), which indicates the average strength of the geostrophic currents (Hermans et al., 2020a). The pattern of the interannual SSH variability is compared with a satellite altimetry sea level anomaly product.

Figure 5 Comparison of modelled SSH with satellite altimetry (mm). (a) Ensemble mean SSH anomaly (2000–2019), (b) ensemble interannual variability (standard deviations of annual mean SSH, averaged over the ensemble, 2000–2019), (c) AVISO MDT (mean SSH), (d) CS3 interannual variability of sea level anomaly.

We use the AVISO CLS18 MDT product (Fig. 5c), estimated from the period 1993–2012 (Mulet et al., 2021), and compare it to the NWSPPE ensemble mean SSH for the present-day ensemble statistics (2000–2019). As AVISO has a lower horizontal resolution (0.25°) and is not able to resolve features with spatial wavelengths of less than ∼ 180 km (Legeais, 2018), we smooth the model by convolving with a uniform filter of 13-by-13 grid boxes (∼ 90 km, Fig. 5a).

We compare the simulated interannual SSH variability 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 (temporally), and calculate the temporal standard deviation. This is compared with the interannual variability from the present-day ensemble statistics (2000–2019) – the square root of ${\mathit{\sigma }}_{\text{int}}^{\mathrm{2}}$ from Eq. (A1) in Sect. A.4.3.

Figure 6 Tide gauge–model sea level anomaly comparison (metres). (a–t) Time series of annual mean SSH from tide gauges (red) and the 12 ensemble members (grey), with the ensemble mean ± 2 ensemble standard deviations shown in black. No account is taken of the differing benchmarks between the model and the tide gauges, so an arbitrary offset is added to align the two datasets to allow for visual comparison of the trends and variability. Note that the differing tide gauge record lengths lead to different ranges on the x and y axes. (u) Tide gauge locations are given in Table 3.

Table 3 Tide gauges used for model comparison.

3.5 Tide gauge data

The NWSPPE was compared with 20 tide gauges around the UK and NWS (Fig. 6, Table 3), which were selected for the length of overlap with the NWSPPE. These were downloaded from the Permanent Service for Mean Sea Level (PSMSL; Holgate et al. (2013); data retrieved 24 November 2022) and processed into annual means. 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 are data from at least 11 months.

The SSH from the nearest (sea) grid box was extracted from the monthly mean data and processed into annual means for each ensemble member and then processed into ensemble means. The tide gauges were not corrected for glacial isostatic adjustment (GIA), and so Smögen (which exhibits substantial GIA) was excluded. For each tide gauge, an offset between the ensemble mean was removed. This was the difference between the NWSPPE ensemble mean temporal mean and tide gauge temporal mean for the common overlap period.

Table 4 Evaluation of transport cross-sections, giving the cross-section name and number, the longitude and latitude of the two ends, the direction that is considered positive, and whether the transport is the positive component or the net transport. The observed estimates (with error estimates where available) and the NWSPPE ensemble mean annual mean (and seasonal range) are given.

Figure 7 Comparison of modelled and observed volume transport through cross-sections. (a–l) The NWSPPE ensemble mean (bold black) seasonal cycle (and ensemble members, coloured thin lines) compared with observed transport estimate (red circle, and horizontal line) and observation error estimate (vertical red line) where available. (m) Map of the NWS giving the location of the cross-sections. Table 4 gives the longitudes and latitudes of the ends of the cross-sections and summarise the data.

3.6 Volume transport estimates

We assess the NWSPPE volume transport through several cross-sections compared to observational estimates from the literature (Fernand et al., 2006; Brown et al., 1999; Svendsen et al., 1991; Turrell et al., 1992; Danielssen et al., 1997; Prandle et al., 1996; Holt et al., 2001), following Lowe et al. (2009). We define these cross-sections in Table 4. Most observational estimates are for the net transport, apart from cross-sections 1, 4, 5 and 12 (Fig. 7m), which are for the positive transport component. Most modelled transport cross-section time series have been processed post hoc and so are limited to net transport, with sections 1, 4, 5 and 6 being the exceptions, where the net as well as positive and negative transport components are available. Section 12 (Shelf Current FS Current) is based on the observed positive (northwestward) transport component and compared with the net modelled transport; therefore, the modelled estimates include a negative component, which would reduce the magnitude of the modelled transport.

4 Evaluation

We evaluate the NWSPPE between 2000 and 2019 for tides, temperature, salinity, sea surface height and transport (Table 1). We show the modelled cotidal charts and compare these to other modelling systems. We compare the modelled SST with the OSTIA analysis (Roberts-Jones et al., 2012), and temperature and salinity with the EN4 profiles dataset (Good et al., 2013) and ICES climate time series (González-Pola et al., 2022). We compare the modelled sea surface height to tide gauge data retrieved from the Permanent Service for Mean Sea Level (PSMSL; Holgate et al., 2013) and two satellite products (Rio et al., 2014; Legeais et al., 2018). We compare modelled transports to observation-based estimates reported in the literature.

Evaluating climate projections is always complicated by the fact that they are not designed to simulate the observed phase of weather and climate variability. For example, while the approximate number and frequency of warm years should agree with observations, their ordering will not. This means the climate must be evaluated rather than the weather. To do this, we compare the statistics of long records of observations with the models – we would expect a level of agreement between the 20-year mean (and standard deviation) of the model and the observations. This works well for SST, where there are satellite records since the 1990s, but other long time series of observations are relatively sparse. Another complication is that reality is a single realisation of what conditions could be expected, given the current state of the climate. If there were 100 worlds with today's greenhouse gas and aerosol conditions, we would expect a range of temperatures, due to unforced variability – different states of the North Atlantic Oscillation, different states of the Atlantic Multidecadal Variability and different weather conditions. With our PPE, we have a range of realisations of the present-day climate, any of which could be the reality. We therefore ask if the observed climatology is consistent with the distribution of modelled climatologies (the NWSPPE). We do this by calculating the observed 20-year mean and the modelled 20-year mean for each ensemble member. We then calculate the normalised bias using Eq. (1) to ask how many (PPE ensemble) standard deviations the observed 20-year mean is from the PPE ensemble mean. If the (absolute) normalised bias is less than 2, the observations are within 2 standard deviations of the ensemble mean, and we consider the PPE to be consistent with the observations.

The approach works well for OSTIA SST and the ICES time series, where there are regular observations every year between 2000–2019. It is complicated by the sparsity of the EN4 data, where there may only be a single observation for a given grid box for a 20-year period. This means the observed weather is not separated from the observed climate – this must be taken into consideration (see Sect. 3.2.1).

4.1 Tidal evaluation

We have evaluated the tides by producing cotidal charts for leading constituents. The tidal harmonic analysis is done online by CO9 for consecutive 20-year periods. We show the ensemble mean for the files created for the 2000–2019 period (Fig. 1). The ensemble mean phase (angle) is converted to its northward and eastward components before averaging. The M₂ component visually agrees well with O'Dea et al. (2012, 2017).

4.2 Sea surface temperature

We find the NWSPPE ensemble mean SST has a high spatial correlation (> 0.95) with the OSTIA analysis across the domain (for all months, seasons and in the annual mean). Figure 2 gives the absolute SST bias (NWSPPE ensemble mean minus OSTIA SST) for the annual mean and the four seasons. When averaged over the domain, the absolute mean bias is less than 0.4 °C for all seasons (and months and in the annual mean). We also assess where OSTIA sits within the distribution of the NWSPPE and indicate where the OSTIA SSTs are more than 2 standard deviations from the ensemble mean with cross-hatching in Fig. 2.

In the annual mean, we find that the OSTIA SSTs of most of the NWS are within the ensemble spread of the NWSPPE; 95 % of the NWS have OSTIA SST biases less than 0.68 °C. Around the edge of the NWS (as delineated by the green contour in Fig. 2), the shelf slope current is slightly warmer than the ensemble, and the Norwegian Trench is slightly cooler. In the wider oceanic parts of the domain, the Norway Basin (63° N, 5° W) is too cold, and a region along the western boundary that extends to 15° W is too warm.

When looking across the seasons, we see that OSTIA SSTs tend to be within the NWSPPE, albeit with less agreement than in the annual mean. While the winter and autumn OSTIA SSTs are largely within the NWSPPE, the summer OSTIA SSTs to the west of the UK are warmer than OSTIA, while the spring OSTIA SSTs in the northern North Sea and around Scotland are cooler than the NWSPPE. Despite the regions of warm summer bias and cool spring bias, overall we consider the NWSPPE SSTs to be broadly consistent with the OSTIA analysis.

4.3 Subsurface temperature and salinity

The EN4 data processing procedure (outlined in Sect. 3.2.1) produces sparsely sampled maps of the NWS, at monthly, seasonal or annual granularity. The seasonal maps are given in the Supplement (SST: Fig. S4; SSS: Fig. S5; NBT: Fig. S6); however, as they are still relatively sparse, they are hard to interpret quantitatively. We have therefore averaged these data across the evaluation regions (in green in Fig. S3) and presented the regional means and spread (2 time the spatial standard deviation) for the larger validation region mask in Fig. 3. For completeness, we also include the regional summary for the Wakelin et al. (2012) region mask in Fig. S3.

This shows that, in most regions and in most times of the year, the EN4 SST observations are within the PPE, although less so in the poorly sampled Norwegian Trench. The NBTs also tend to be within the PPE and tend to agree with the SST in most regions and months. There is more disagreement with the EN4 SSS, with much greater spatial variability within the regions (reflecting the greater spatial variability in Fig. S6). However, when averaging over the regions, the EN4 SSS is typically within the PPE.

4.4 In situ temperature and salinity time series

We now evaluate the NWSPPE against the ICES climate time series. All ICES time series have absolute temperature biases less than 1.25 °C (NWS time series < 0.85 °C) and absolute salinity biases of less than 0.75 (Fig. 4, Table S1). Most of the NWS ICES temperature time series have a good overlap between the model and the observations mean values, with NWS in situ observations being within the NWSPPE. The Fair Isle Current is modelled too cold, and outside the NWSPPE, although the absolute temperature bias is −0.83 °C. All off-shelf ICES temperature time series are too cold and are outside the NWSPPE, reflecting the region of cold bias in the Norway Basin shown in the OSTIA analysis. Conversely, all the off-shelf ICES salinity time series are within the NWSPPE, with very low absolute biases (< 0.15), suggesting a good agreement. The NWS ICES salinity time series are typically too salty, which may, in part, relate to the coastal locations of many of the sites, where there are typically large horizontal salinity gradients. The Fair Isle Current is 0.42 too salty and is outside the NWSPPE (with a large relative bias) with little overlap of the interannual variability between the observations and model.

We can also use the ICES time series to assess the interannual variability. On the shelf, most locations have similar interannual variability between the observations and the model, with the observations sitting within the spread of the ensemble (having an absolute Relative Interannual Variability Bias less than 2). The Plymouth WCO E1 temperature and Helgoland Roads salinity are the exceptions, with the observations being much more variable than the model. Off the shelf, the observed salinity interannual variability is well represented in the NWSPPE, while the Faroe Shetland NAW temperature interannual variability is modelled as too variable compared with the observations.

Overall, the NWSPPE is in relatively good agreement with the ICES time series, giving further confidence to the NWSPPE simulations.

4.5 Sea surface height

COx has a non-linear free surface, and so it simulates the dynamic response of the sea level to the local dynamics of the model. This illuminates important aspects of the model behaviour. It is also a component of sea level but should not be used directly as a set of sea level projections (see Sect. A.5).

First, we compare the NWSPPE and the AVISO CLS2018 satellite MDT (Fig. 5). There is a good agreement in the overall pattern, with a large-scale NW/SE gradient and a similar range of values. On the NWS, the highest levels are in the German Bight, with lower values in the northern North Sea, north of the Dooley Current.

We then compare the simulated interannual SSH variability to that of the C3S SLA product (Legeais et al., 2018). There is good spatial agreement between the PPE and C3S SLA product. In the open ocean, there is greater sea level variability in deeper regions (e.g. the Rockall trough and the Icelandic basin). On the shelf, the greatest variability in both the PPE and the C3S SLA product is in the German Bight, and the lowest variability is in the Celtic Sea and to the west of Ireland and adjacent the shelf break; this is more pronounced in the C3S SLA product.

There is good agreement in the trend and interannual variability of the tide gauges and the PPE (Fig. 6). The interannual variability of the PPE ensemble mean is averaged out, but the tide gauge variability tends to agree with the ensemble spread (ensemble mean ± 2 standard deviations). In most locations (apart from Dieppe and Bergen), there is a good agreement in the trends. The tide gauge at Dieppe was re-installed in 2009, and the levels appear to have changed (https://psmsl.org/data/obtaining/stations/474.php, last access: 20 May 2024), with the later values in closer agreement with nearby tide gauges (La Havre and Boulogne). The tide gauge at Bergen may be affected by GIA.

4.6 Volume transport through cross-sections

We now evaluate the NWSPPE residual circulation. While it is difficult to directly measure residual currents (as they are typically much smaller than the oscillatory tidal currents), there are several observation-based estimates of volume transport through cross-sections across the NWS (Fig. 7m, Table 4) that may be used for model evaluation (Fernand et al., 2006; Brown et al., 1999; Svendsen et al., 1991; Turrell et al., 1992; Danielssen et al., 1997; Prandle et al., 1996; Holt et al., 2001). We use the same 12 cross-sections as used by Lowe et al. (2009) and Tinker et al. (2015). Here we report the modelled transport as an ensemble annual mean and the seasonal range of the ensemble mean (minimum–maximum)

Overall, there is a relatively good agreement between the model and the observations, with all transports being of the correct direction and with the observations generally being within the modelled seasonal cycle. Furthermore, the transports through most cross-sections are similar or show an improvement when compared with the UKCP09 (Lowe et al., 2009) and Minerva (Tinker et al., 2015) projections.

The North Sea inflow through the Fair Isle (6) is modelled to be too strong (0.32 Sv, 0.20–0.46 Sv) compared with observations (0.2 Sv), while the outflow through the Norwegian Trench (5; 1.56 Sv, 1.31–1.81 Sv) is modelled as being too weak (1.8 Sv), although the observations are overlapping with the modelled seasonal cycle. The North Sea inflow between Shetland and the Norwegian Trench (4: 0.61 Sv, 0.51–0.71 Sv) and through the Dover Strait (10: 0.07 Sv, 0.03–0.12 Sv) are in good agreement with the observations (0.60 Sv and 0.1 Sv, respectively).

There is a large inflow and outflow from the North Sea in the Skagerrak, which are nearly equal, with a very small residual net outflow. This modelled inflow (1; 1.21 Sv, 0.94–1.41 Sv) is very close to the observations (1.0 Sv, 0.5–1.5 Sv). The Dooley Current (3; 0.25 Sv, 0.23–0.30 Sv) and the eastward flow to the north of the Dogger Bank (2; 0.08 Sv, 0.03–0.11 Sv) are in good agreement with the observations (0.25 Sv, 0.12–0.38 Sv and 0.05 Sv, respectively). The St George outflow (9: 0.11 Sv, 0.02–0.20 Sv) is in good agreement with the observations (0.18 Sv). While the Irish shelf (8: 0.12 Sv, −0.04–0.20 Sv) and the Hebrides shelf (7: 0.43 Sv, 0.33–0.58 Sv) are modelled as being too weak (0.25 Sv) and too strong, respectively (0.25 Sv), although where cross-sections are not across channels, the transport is sensitive to the exact width of the cross-section. The Shelf Current FS Current (12: 2.68 Sv, 1.08–3.89 Sv) is modelled as being too weak (5.50 Sv, 4.00–7.00 Sv); however, this could be due to the comparison of the modelled net transport and the observed northwestward transport. Overall, the transport comparisons suggest the NWSPPE simulates the configuration of the NWS residual transport well, and it generally captures the magnitude of the transport pathways.

5 Results and applications

We now explore some of the results of the NWSPPE climate projections – we refer the reader to Tinker et al. (2020) for a description of the results of the present-day control simulation. We first describe the changes to the residual circulation, before considering the spatial patterns of the climate projections and then turning to the temporal evolution with regional mean time series.

5.1 Changes in circulation

We assess the NWS residual circulation to give a background to the NWS changes. We consider how the barotropic residual currents change across the 21st century (2000–2019 and 2079–2098), as well as transport through pre-defined cross-sections (Figs. S7–S14, Table S2).

Figure 8 Barotropic residual current statistics. The early 21st century (2000–2019), late 21st century (2079–2098) and difference of the ensemble mean current magnitude (a, e and i, respectively); variance, as measured by the area of the 2.45 standard deviation current ellipse (b, f and j); the ensemble OVL (the overlap of the ensemble, as a measure of how similar the ensemble members are (c, g and k)); and the early and late 21st century current configuration (d, h; colours show the direction of the currents where the mean is more than 1 standard deviation); l for each ensemble member, the OVL is computed between the early and late 21st century – here we present the ensemble mean of these OVLs.

5.1.1 Mean residual circulation

There is a general weakening of most of the NWS circulation over the 21st century (Fig. 8i), but the general configuration of the residual circulation is relatively consistent across the ensemble and the 21st century (Fig. 8d and h). A small region that extends from the south to the west of Ireland is an exception to both points, which we consider later. There is no evidence of a change in the large-scale configuration of the northern North Sea circulation as seen by Tinker et al. (2016) and analysed by Holt et al. (2018). However, to the west of Ireland (52° N, 13° W), we do see a slight change in configuration related to a localised change in the slope current.

There is a general increase in the strength of the shelf slope current south of 54° N (Figs. 8i and S7) and a decrease between 54° N and 60° N at 6° W (Figs. 8i and S8). To the north of the North Sea (northwest of 60° N, 6° W), there is a substantial strengthening of the slope current (Fig. S9). The portion of the shelf slope current that follows the isobaths turning into the Norwegian Trench slightly increases from ∼ 0.8 to 1.2 Sv over the 21st century (Fig. S10b); however, as the strength of the slope current increases, this represents a relative decrease from 35 % to 10 % (Fig. S10b cf. Fig. S10a). Similarly, the Norwegian Trench outflow is similar (Fig. S11c) but represents a much smaller proportion of the total downstream net transport (40 % to 20 %, Fig. S11c, cf. Fig. S10b). This is consistent with the finding of Holt et al. (2018), who showed that an increase in haline stratification in this region increases the Rossby radius, and so this makes it harder for the current to follow the bathymetry. When looking across the ensemble, all members show a strengthening of the most northern sections of these currents (north of 61° N, 3° W), but there is a general decrease of the currents further south, in the main part of the Norwegian Trench (Fig. S15g and h).

There is an overall weakening of the northern North Sea inflow (NNSI), the Norwegian Trench outflow (Fig. S15h and i) and the wider North Sea circulation. The current though Pentland Firth (58°45^′ N, 3°6^′ W; Fig. S12a) and the Fair Isle Current (59°30^′ N, 1°45^′ W, 0.3–0.2 Sv; Fig. S12b) reduce by about 15 % and 40 %, respectively. The inflow to the west of Shetland increases slightly (14 %, Fig. S12c), but this may reflect an increasing flow into the Norwegian Trench. The Dooley Current also shows a 20 % reduction in the ensemble mean (from 0.256 to 0.202 Sv, Figs. S14a and S15j–l). Further south, the North Sea inflow though the English Channel decreases by 40 % (from 0.07 to 0.04 Sv, Fig. S13a), and this is reflected in the lower downstream currents though the Southern Bight (Fig. S13b) towards the Skagerrak. The Skagerrak recirculation weakens from 1.2 to 1.0 Sv (Fig. S14bc), although this reduces to 0.7 Sv for one ensemble member.

There is a tendency for an increase in the residual current magnitude to the southwest of the UK and Ireland, which leads to a slight change in circulation configuration: 80 % of this region shows an increase in the residual current magnitudes, whereas 60 % of the rest of the NWS shows a weakening. In the present day, there is a clockwise coastal current around Ireland and a northward slope current that follows the shelf break from the Goban Spur (49°30^′ N, 12° W), around the Porcupine Seabight (50.5° N, 12° W) and Porcupine Bank (52° N, 14° W), before following the shelf break north around Ireland. Between the coastal current and the shelf slope current, there is a weak, disorganised southward countercurrent. In the future, the coastal current weakens and the slope current strengthens, which helps strengthen this part of the slope current. The southward countercurrent strengthens substantially by the end of the 21st century, and it flows southeastward across the Celtic Sea, increasing the mean current magnitude and reducing its (relative) variability.

5.1.2 Residual current variability

We now consider the residual current uncertainty ellipse approach of Tinker et al. (2022). For each grid box, we can fit ellipses around the bivariate normally distributed residual currents to encapsulate 95 % of the data – we can then use the properties of these ellipses, and how they vary over time and the ensemble, to describe the residual current similarity. To aid the reader in the interpretation of residual current uncertainty ellipses, we have included examples focusing on the regions with the greatest changes (Figs. S15 and S16). Furthermore, Tinker et al. (2022) include explanatory figures in their appendices (reproduced in our Fig. S2) and released a Python toolbox to undertake the analysis.

Over most of the shelf, there are only small percentage changes in variability of the residual currents (e.g. Figs. S15 and S16), as given by the uncertainty ellipse area (Fig. 8b, f, and j). While there are no NWS regions with a large reduction in variability, there are substantial regions with large variability increases (e.g. Fig. S16d, e and g) – these include the NNSI (to the west of Ireland), the slope current and the Norwegian Trench. Where this increase in variability is coupled by a decrease in the mean current magnitude, we see a large decrease in the significance of the residual circulation (such as in the NNSI, the Dooley Current and the Norwegian Trench). Where the increase in variability is matched by an increase in the mean current magnitude, the significance of the current tends to increase (i.e. the currents to the west of Ireland), but this is not always the case (e.g. the slope current to the west of Shetland).

The residual current ellipses allow us to compare the variability of the residual current to its mean and quantify how many standard deviations the mean is from zero; when this is greater than 2.45, the mean current does not cross zero more than 95 % of the time. As we have shown that the residual circulation magnitude is typically reducing and that the variability is increasing, the number of standard deviations is also generally decreasing, and this is reflected in the amount of NWS residual circulation that is significant (Fig. 8d and h), which we use to show the NWS residual circulation configuration.

We can also use the residual current ellipses to quantify similarity between two periods, between two ensemble members or across an ensemble. As each residual current ellipse represents a bivariate normal distribution, which integrates to 1, we can integrate the volume under two or more bivariate normal distributions to show how similar they are (illustrated in Fig. S2i), where 1 is perfect agreement, which decreases towards zero as the mean, variance or covariance of the distributions differs. This allows us to quantify the residual circulation ensemble spread in the early and late 21st century and how different the residual circulation is between the early and late 21st century.

In addition to the general reduction in the ensemble mean residual current magnitude (Fig. 8i) and increase in residual current variability, there is also an increase in the ensemble spread and diversity. This is reflected in the decrease in ellipse ensemble overlap (Fig. 8c, g, and k), which is a measure of how similar the pointwise residual current distributions are across the ensemble. When this is high for a given point (near the maximum value of 1), the residual current distribution (in terms of U and V mean and variance, as well as their covariance) are similar across the ensemble. In the early 21st century (Fig. 8c), over half of the NWS has and ens_OVL > 0.65, whereas by the late 21st century there has been a general reduction across most of the NWS with ens_OVL < 0.53 (Fig. 8g). But there is a particular reduction in the northern North Sea, Norwegian Trench, and to the south and west of Ireland (Fig. 8k), where ens_OVL < 0.3. Conversely, there is a slight increase in ens_OVL in the Irish Sea, suggesting a convergence in the ensemble behaviour.

ens_OVL quantifies the spread of the ensemble in the early and late century (Fig. 8c and g), as well as how much it has changed (Fig. 8k). It is possible for ensemble spread to remain the same (ensemble OVL does not change over the 21st century) but with the future being nonetheless very different from the present day. To quantify how much the circulation changes over the 21st century in each ensemble member, we calculate the ellipse overlap between the early and late 21st century, and we report the ensemble mean in Fig. 8l. This shows that over most of the NWS there is little change in the circulation, whereas the northern North Sea, Norwegian Trench and to the west of Ireland there are regions of substantial change. While this is clear in the ensemble mean, some ensemble members show little change in their circulation, and others are showing a significant change reflecting changes in their current magnitudes and variances.

Figure 9 The SST ensemble mean change over the 21st century (2079–2098 compared with 2000–2019), for the (a) annual mean, (b) winter (December–February), (c) spring (March–May), (d) summer (June–August) and (e) autumn (September–November). Black contours relate to the colour bar ticks, and the zero contour is a thicker black line. The coastline and NWS region are delineated with a grey line. For a full comparison between the early and late century, as well as the components of variance (ensemble variance and interannual variability), see Fig. S17.

Table 5 Projected regional mean SST changes between 2000–2019 and 2079–2098. Ensemble mean changes are given with ± 2 (ensemble) standard deviations.

We have produced some exemplar uncertainty ellipses (for example locations) to help visualise some of these changes. The present-day ellipses tend to be very similar across ensembles, which is reflected by ens_OVL > 0.6 in most locations. The Irish Sea (e.g. Fig. S16a and b) and the outermost point to the west of Ireland (e.g. Fig. S16c) are the exceptions, which are captured in Fig. 8c, where the ensemble tends to converge. In the northern North Sea (Fig. S15) and to the south and west of Ireland (Fig. S16), the uncertainty ellipses are larger in the future, and there is much less agreement across the ensemble, with less overlap, captured by a much lower future ens_OVL (e.g. Fig. S15g). Other regions have a similar level of ensemble agreement (e.g. Fig. S15d) and so less change in the ens_OVL over the 21st century. Some locations show relatively little change in the ensemble spread (similar values of the ens_OVL over the 21st century) but a large difference in the residual currents (e.g. Fig. S15e), reflecting little overlap between the present day and the future, which is captured by a low OVL.

5.2 Projected changes at the end of the 21st century

The NWS SST (Figs. 9 and S17; Table 5) rises by 3.11 °C (± 0.98 °C) by the end of the century (2079–2098 compared with 2000–2019), when averaged across the shelf (ensemble mean ± 2 ensemble standard deviations). There is a seasonal cycle in the warming, with greater warming in the summer (3.57 °C ± 1.09 °C) and autumn (3.73 °C ± 1.07 °C) than in the winter (2.72 °C ± 0.97 °C) and spring (2.43 °C ± 1.01 °C). This is consistent with the wider UKCP18 projections (Murphy et al., 2018), which exhibit greater warming in the summer than the winter, which increases the amplitude of the temperature seasonal cycles, and is consistent with their key findings of projected hotter drier summers and warmer wetter winters (Lowe et al., 2019). There is a region of reduced warming to the NW of the NWS (e.g. 60° N, 15° W) which is more visible in winter (it may be obscured by summer stratification) and may influence the NWS. A likely cause of this is the North Atlantic warming hole (Menary and Wood, 2018) and a reduction in the Atlantic Meridional Overturning Circulation (AMOC), transporting less heat from the tropical Atlantic. This will be considered later in the Discussion section.

On the NWS, the southern and central North Sea show considerable warming, particularly in summer and autumn, with values typically greater than 3.5 °C warming. The greatest NWS warming is in the Dover Strait in autumn, where the ensemble mean SST warms by up to 5 °C. There is an adjacent region of reduced warming in the centre of the Southern Bight (∼ 52° N, 3° E) in all seasons, relating to an apparent weakening of the warm plume of water flowing from the English Channel into the southern North Sea (Fig. S13a).

Figure 10 The NBT ensemble mean change over the 21st century (2079–2098 compared with 2000–2019), for the (a) annual mean, (b) winter (December–February), (c) spring (March–May), (d) summer (June–August) and (e) autumn (September–November). Black contours relate to the colour-bar ticks, and the zero contour is a thicker black line. The coastline and NWS region are delineated with a grey line. For a full comparison between the early and late century, as well as the components of variance (ensemble variance and interannual variability), see Fig. S18.

Table 6 Projected regional mean NBT changes between 2000–2019 and 2079–2098. Ensemble mean changes are given with ± 2 (ensemble) standard deviations.

The near-bottom temperatures exhibit a more modest rise than the SST (Figs. 10 and S18; Table 6), with an annual mean NWS NBT warming of 2.49 °C (± 0.94 °C). This largely reflects a reduced warming under stratified regions. For the SST, we see a greater warming in summer than winter. As the spring initialisation of stratification isolates bottom water from the atmosphere, NBTs remain close to the winter and spring temperatures under stratification. Therefore, the weaker winter SST warming will reduce the summer NBT warming in stratified regions. This is illustrated in the difference in the southern and northern North Sea NBT summer warming: 3.55 °C (± 1.00 °C) and 2.07 °C (± 1.00 °C), respectively. The southern North Sea summer NBT rises to a similar level to that of the SSTs in the same region and season (3.78 °C (± 1.04 °C)), while the northern North Sea summer NBT warming is closer to the winter SSTs (2.71 °C (± 0.99 °C)) than the summer SSTs (3.61 °C (± 1.20 °C)) in the same region.

It is also interesting to note that the northern North Sea winter (and autumn) NBTs exhibit greater warming than the summer. The northern North Sea is the main place where the North Sea connects to the North Atlantic, and so it is most affected by the reduced North Atlantic warming. In the winter, the northern North Sea may also be warmed by the atmosphere, while when stratified this would not be possible. This would tend to lead to a greater northern North Sea NBT rise in winter than in summer, as simulated.

Figure 11 The surface–bed temperature difference (DFT) ensemble mean change over the 21st century (2079–2098 compared with 2000–2019) for the (a) annual mean, (b) winter (December–February), (c) spring (March–May), (d) summer (June–August) and (e) autumn (September–November). Black contours relate to the colour-bar ticks, and the zero contour is a thicker black line. The coastline and NWS region are delineated with a grey line. For a full comparison between the early and late century as well as the components of variance (ensemble variance and interannual variability), see Fig. S19.

Table 7 Projected regional mean DFT (the difference between SST and NBT) changes between 2000–2019 and 2079–2098. Ensemble mean changes are given with ± 2 (ensemble) standard deviations.

The fact that the stratified SSTs tend to rise more than the NBTs, leads to an increase in the difference between the surface and bed temperatures (DFT, Figs. 11 and S19; Table 7), which reflects an increase in the stratification magnitude. Stratification can be quantified by the amount of energy needed to mix the stratified water column, which is formalised with the potential energy anomaly (PEA), defined in Eq. (2). There is a general strengthening in the magnitude of summer stratification, although there is little change in its spatial extent (Fig. S21, Table S4). PEA can be broken down into thermal and haline components (PEAT and PEAS, Table S5 and S6, respectively). Most of the increase in PEA is associated with an increase in thermal stratification (PEAT), although in places this is modified by changes in the haline stratification (PEAS). Around the shelf edge between the Atlantic and the NWS, there is a strengthening in the haline stratification throughout the year. This balances a co-located band of reduced thermal stratification (PEAT) in the winter and spring. This haline stratification extends into the western half of the Celtic Sea, leading to a substantial increase in local winter PEA. Finally, there is a reduction in haline stratification adjacent to the Norwegian Trench, which at its peak in summer balances any increase in PEAT, leading to no increase in PEA at ∼ 60° N, 2.5° E.

The changes we see to the NWS temperatures also reflect a change in the seasonal cycle. The SST seasonal cycle subtly changes between the early and late 21st century (when averaged over the NWS), with a slight preferential warming between late June and early December, which leads to the greatest surface warming being in autumn. This is consistent with the broader behaviour of the PPE and (a subset of) CMIP6 models (Cotterill et al., 2023), which shows an atmospheric circulation change that increases the frequency of drier summer-type regimes and a decrease in stormy winter types in autumn. There is little change in the NBT seasonal cycle between the early and late 21st century (when averaged over the NWS). As the NWS is still mixed in the winter throughout the 21st century, the change in the seasonal cycle of DFT (the difference between the surface and bed temperatures) is dramatic. From January to May, DFT (averaged over the NWS) does not change over the 21st century, while at its peak (in mid-August) NWS mean DFT increases by 1.5 °C. These seasonal maximum increase in NWS DFT over the 21st century is approximately a month later than the present-day NWS DFT seasonal peak. Therefore, there is little difference to the DFT in the early summer and a much greater change in the autumn, which effectively extends the seasonal stratification.

Figure 12 The ensemble mean duration, initialisation and breakdown of stratification in day of the year (first, second and third columns, respectively) for the early 21st century (2000–2019), late 21st century (2079–2098) and how it changes between the two periods (first, second and third rows, respectively).

Table 8 Projected regional mean SSS changes between 2000–2019 and 2079–2098. Ensemble mean changes are given with ± 2 (ensemble) standard deviations. Salinities are given on the (unitless) practical salinity scale.

We assess the timing of onset and breakdown of the stratification, as well as its duration, using the method described in Sect. 2.9. We find that average NWS onset is typically 4.16 (± 2.31) d earlier when averaged over the ensemble (± 1 ensemble standard deviation), and breakdown is 18.91 (± 3.50) d later, leading to an increase in stratification of 23.07 (± 5.06) d. This is a greater increase in stratification duration than reported by Holt et al. (2010), who found that the onset of stratification was 1.29 d earlier and the breakdown was 3 d later, giving a much smaller increase in duration of 6.36 d when averaged over the NWS (Fig. 12).

Figure 13 The SSS ensemble mean change over the 21st century (2079–2098 compared with 2000–2019) for the (a) annual mean, (b) winter (December–February), (c) spring (March–May), (d) summer (June–August) and (e) autumn (September–November). Black contours relate to the colour-bar ticks, and the zero contour is a thicker black line. The coastline and NWS region are delineated with a grey line. For a full comparison between the early and late century, as well as the components of variance (ensemble variance and interannual variability), see Fig. S20.

There is a substantial freshening across the domain, which is greatest in the open ocean to the west of the domain (Figs. 13 and S20; Table 8). The NWS also shows substantial freshening, with −1.01 (± 0.93) when averaged over the shelf, and little spatial or seasonal variations. There also is a substantial increase in the ensemble variance and a much smaller increase in interannual variability. This increase in salinity appears to be oceanic in origin and predominantly entering the NWS via the NNSI. Regions that are less connected to the NNSI show a slightly reduced freshening. For example, there is slightly less freshening in the Celtic Sea, English Channel and southern North Sea, which in part may reflect the use of climatological river forcings but also the more southernly origin of their Atlantic water. There is also a reduced freshening in the Norwegian Trench, which is likely related to the use of climatological Baltic lateral boundary conditions (LBCs).

When looking at surface salinity of the wider North Atlantic in the HadGEM3-GC3.05 PPE, we start to see a meridional freshening dipole, with a general freshening north of ∼ 40° N (20° N in the Canary Current) and with an increase in salinity to the south (Fig. S23). This pattern strengthens throughout the 21st century and across the ensemble. The greatest freshening is from the Gulf of St Lawrence and the Labrador Sea, as well as along the pathway of the Gulf Stream and North Atlantic Drift, towards the NWS. There are regions of increasing salinity north of the dipole, including Hudson Bay and in the Norway Basin (63° N, 5° W). The same pattern and temporal evolution occurs in all ensemble members, although the timing and strength of the pattern varies. The substantial North Atlantic freshening across the ensemble is highly correlated with the NWS freshening (Fig. S24a) and is of a similar strength (Fig. S24b), suggesting the NWS response is driven by the Atlantic and that it also explains the increase in ensemble spread.

Figure 14 Temporal evolution of the regional mean annual SST. Each region is given in a separate panel (a–n), with panel (o) showing the NWS region mask. Each ensemble member is filtered with a 5-year low-pass filter (shown faintly), with the ensemble mean and ± 2 standard deviations shown in red (both based on the filtered ensemble members). The Copernicus Marine RAN reanalysis regional mean SST (also filtered) is shown in black for comparison. The y-axis limits change with each region.

Figure 15 Temporal evolution of the regional mean annual SSS. Each region is given in a separate panel (a–n), with panel (o) showing the NWS region mask. Each ensemble member is filtered with a 5-year low-pass filter (shown faintly), with the ensemble mean, and ± 2 standard deviations shown in red (both based on the filtered ensemble members). The Copernicus Marine RAN reanalysis regional mean SST (also filtered) is shown in black for comparison. The y-axis limits change with each region.

5.3 Temporal and ensemble evolution over the 21st century

We now assess the regional mean time series from the NWSPPE to illuminate the NWS's temporal evolution. We also compare to the Copernicus Marine reanalysis (RAN) to provide near-present-day information for context (Figs. 14 and 15). We note that our use of climatological Baltic boundary conditions may explain the behaviour and the reduced interannual variability in the Skagerrak/Kattegat and Norwegian Trench regions.

The SST regional mean time series are in relatively good agreement with the absolute values of the RAN (Fig. 14). There is a near-linear trend in the ensemble mean SST. Each individual ensemble member also has near-linear trends (typically r > 0.9) but have much larger interannual variability, which is averaged out in the ensemble mean (Fig. S25). While the ensemble mean interannual variability appears much lower than that of RAN, there is relatively good agreement when comparing to individual ensemble members (Fig. S25), particularly on the shelf. The RAN time series appear to have a different trend to the ensemble, but this is consistent with the scale of the variability of individual ensemble members (not shown). The SST ensemble spread remains relatively constant in some regions and slightly widens in other regions. The SST interannual variability is also relatively constant throughout the projection (Fig. S25). The annual mean NBT data show similar behaviour to that of the SST and so are not included in this paper.

The SSS regional mean time series tend to be saltier than the RAN reanalysis, which lies outside the ensemble spread in several regions. This is particularly true in the Skagerrak and (downstream) in the Norwegian Trench where the NWSPPE is much fresher than the RAN. This difference partly reflects the differences in the Baltic lateral boundary conditions between the RAN and the NWSPPE, as well as the complexity of modelling this region, especially when using Baltic and Atlantic lateral boundary conditions from a different ocean model (as the RAN and NWSPPE both do). Moreover, the RAN surface salinity to too salty in the Norwegian Trench (the Skagerrak is outside the RAN evaluation region) (Renshaw et al., 2019).

There is good agreement in SSS interannual variability with the reanalysis in the North Sea and English Channel, while the PPE has greater interannual variability in other shelf regions. There is also a substantial increase in SSS interannual variability between 1990–2019 and 2069–2098 in almost all regions of the NWS (Fig. S26) – absent in SST (Fig. S25)

The NWS SSS behaves very differently to that of SST (Fig. 14). Rather than the near-linear evolution seen in SST (Fig. S27), SSS has an increasing freshening and divergence of the ensemble. For most ensemble members and shelf regions, the correlation tends to improve when fitting different straight lines for the beginning and end parts of the time series (Fig. S28). While most ensemble member have similar initial SSS gradients, there is a greater range of gradients in the second half of the century, and a range of timings of the change in slope – both of which contribute to the divergence of the ensemble. This is consistent with the arrival and strength of the freshening signal from the northeast Atlantic.

6 Discussion

We have downscaled a PPE of HadGEM3-GC3.05 to give a set of climate projections for the NWS. The PPE was developed for the UKCP18 climate projections to explore the range of an important source of uncertainty in the climate projections and was run under the RCP8.5 scenario. We have downscaled 12 members of this PPE with the shelf sea model NEMO CO9 (NWSPPE) from 1990–2098 (plus a discarded 10-year spin-up period). We have driven CO9 with atmospheric fluxes and ocean lateral boundary conditions interpolated directly from HadGEM3-GC3.05 (without any bias correction). We have used a climatology for the river forcings and for the exchange with the Baltic Sea. We have undertaken extensive model evaluation (for SST, subsurface T and S, SSH, and tides) and explored the first-order change in the NWS. These simulations represent the state of the art for NWS marine projections.

Figure 16 Comparison of (a) the SST anomaly and (b) SSS anomaly of the current NWSPPE (red) with the Minerva projections (blue) for the NWS. Each ensemble member is filtered, and the 2000–2019 mean is removed to give the anomaly. Note the current ensemble is run under RCP8.5, while Minerva was run under SRES A1B.

This set of projections is closely related to, and compatible with, the 200-year present-day control simulation of Tinker et al. (2020) (PDCtrl). We have combined both datasets into a single data release, which is also described here.

Our NWSPPE is a direct update to the Minerva projections released in 2016. The same general approach is taken, downscaling the PPE run for UKCP, with a regional shelf sea model, to give a similar set of physical projections. Minerva downscaled the UKCP09 PPE which was based on HadCM3 (Gordon et al., 2000; Pope et al., 2000). HadCM3 was a relatively old model at the time (being a CMIP3-generation model) but was used as its computational efficiency allowed it to be run many times in the ensemble. Here, we use the latest Met Office Hadley Centre CMIP6-generation model HadGEM3-GC3.05. Minerva was based on POLCOMS (Holt and James, 2001; Holt et al., 2001), a 12 km resolution shelf sea model which had already been replaced at the Met Office in favour of the NEMO Coastal Ocean model. Here we used NEMO version 4.04, an even more recent version than is used operationally at the Met Office. As the HadCM3 atmosphere has a horizonal resolution of 2.5° latitude by 3.75° longitude, a regional atmosphere climate model (HadRM3, ∼ 25 km horizontal resolution; Jones et al., 2004) was used to downscale the HadCM3 atmosphere to provide atmospheric forcings to POLCOMS. As the HadRM3 model domain did not cover the full POLCOMS model, the southwest corner had to be cut off. HadGEM3-GC3.05 is of a much higher resolution (∼ 60 km) than HadCM3, and so the atmosphere did not need to be downscaled to produce the atmospheric forcing for NEMO.

Our projections are run under the RCP8.5 scenario, while the Minerva projections were based on the SRES A1B scenario (Nakicenovic et al., 2000). There is a greater warming in RCP8.5 than in SRES A1B, with the UK region projected to warm by 3.9 °C in 2080–2099 relative to 1981–2000 under RCP8.5 compared with 2.7 °C under SRES A1B (Lowe et al., 2019) (the 50th percentile of the UKCP18 terrestrial probabilistic projections). Furthermore, all HadGEM3-GC3.05 PPE members have a relatively high climate sensitivity (Rostron et al., 2020), similar to HadGEM3-GC3.1 and greater than that of HadCM3 (Andrews et al., 2019), which was used in the Minerva projections. This greater warming is reflected in our new projections. We project an SST rise of 3.11 °C (± 2σ = 0.98 °C) and an SSS freshening of −1.01 (± 2σ = 0.93) for 2079–2098 relative to 2000–2019. Tinker et al. (2016) projected shelf and annual mean SST rise of 2.90 °C (± 2σ = 0.82 °C) and an SSS freshening of −0.41 (± 2σ = 0.47) for 2069–2098 relative to 1960–1989. It is difficult to directly compare these changes due to the changes in baseline; however, we note that we project a greater change over a shorter period.

When we plot the regional mean time series for both sets of projections, as anomalies relative to a common baseline, the difference is more dramatic (Fig. 16). While both Minerva and the NWSPPE show a near-linear SST rise and similar ensemble divergence, there is a notably greater rate of warming in our new projections. When looking at SSS, the difference is striking. While the Minerva projections showed a substantial reduction in salinity, we show a much greater freshening. In the Minerva projections, the ensemble mean salinity is mainly driven by a step-like change in a few ensemble members with particularly high climate sensitivity. Almost all our ensemble members show a substantial increase in the rate for freshening (Fig. S10), driving a much greater ensemble mean change. Figure 13 suggests the salinity changes are driven by changes in the salinity in the open ocean adjacent to the NWS.

Over most of the NWS, we show a similar pattern of SST change to Minerva (Tinker et al., 2016) with the southern North Sea showing the greatest warming, and the SST fingerprint of a reducing shelf break current to the north and west of Scotland (Fig. 9). When we look across the whole model domain, we see important differences in the SST change spatial pattern (compared with Minerva). The most prominent difference is the area of near-zero winter (and spring) warming in the open ocean to the northwest of the domain in NWSPPE, which is absent in Tinker et al. (2016). This is consistent with a slowdown in the AMOC (Table S9) and the associated North Atlantic warming hole (Drijfhout et al., 2012). Figure S29 (maps of the correlation between AMOC slowdown and NWS change) show a concurrent region of highly correlated SSTs to the south of Iceland and a pattern of highly correlated NBTs in the deep waters to the west of the NWS. These correlations reflect the fact that the ensemble members with the weakest reduction in AMOC have the greatest warming locally and also that the ensemble members with the strongest reduction in AMOC actually have local cooling. The tongue of highly correlated SSTs in the Icelandic Basin is a branch of the North Atlantic Current (Perez et al., 2018), while the pattern of the highly correlated NBTs reflects the complex deep water pathways and bathymetry of this region. There is a competition between the warming associated with a rising global mean temperature and a cooling associated with a reduced AMOC (Drijfhout, 2015); over most of the domain the warming dominates, whereas in this small region the cooling dominates, resulting in this off-shelf region of negligible projected SST warming.

We see a substantial lengthening of the seasonal stratification, which starts 4 d earlier and finishes 19 d later. This is a much greater lengthening than seen by Holt et al. (2010). Stratification is an important feature of the NWS oceanography and has a particularly important impact on the NWS ecosystem. The formation of stratification in spring confines phytoplankton to the sun-lit surface mixed layer of the NWS, providing it with enough light to initialise the spring phytoplankton bloom, which drives the NWS ecosystem. As the summer progresses, the nutrients in the upper layer are used up, and the autumnal breakdown of stratification allows for mixing of subsurface nutrients, leading to a secondary autumn bloom. Most phytoplankton growth happens in the spring bloom, especially diatoms, which dominate in the seasonally stratified zone. There is a much smaller dinoflagellate bloom that happens in late summer/autumn (Van Leeuwen et al., 2015) – our longer projected stratification may suppress mixing and so reduce the nutrient supply to these phytoplankton and restrict their growth; alternatively, it might keep them in the sunlit zone for longer and increase growth, as long as there is a nutrient supply, e.g. from short-lived mixing events. The summer bloom in the stratified North Sea has been observed to be increasing in duration, by extending further into autumn (Silva et al., 2021), although the exact physical mechanism was not identified. Any change in bloom duration or strength would affect primary production and so have a knock-on effect on the rest of the ecosystem. The importance of the spring bloom relative to the autumn bloom may mean that slightly earlier initialisation may be more impactful than substantially delayed autumnal breakdown, as there is potentially more growth there.

We have shown a general weakening of the residual circulation of the NWS, over the 21st century, although its large-scale configuration remains consistent. As the residual currents are important for the transport of heat, salt and matter, their weakening can lead to changes in the mean state. For example, we see a region of reduced warming in the Southern Bight (e.g. Fig. 9). There is eastward transport bringing warm water from the English Channel into the Southern Bight of the North Sea, which leads to a plume of warmer, saltier water (e.g. Tinker et al., 2022). This eastward volume transport decreases into the future (Fig. S13a), as does the associated heat transport. This weakens the Southern Bight warm plume and so reduces the localised warming. Tinker et al. (2016) found a substantial change in the circulation of the northern North Sea and three of their highest climate sensitivity ensemble members, with a reversal in the flow. Holt et al. (2018) undertook an extensive analytical study to replicate this and to understand the mechanisms behind the changes. They reported a change in the North Sea circulation with the shelf slope current bypassing the Norwegian Trench and so reducing the exchange of the North Sea with the North Atlantic, leading to a more estuarine-like circulation for the North Sea. One of the mechanisms they found was an increase in the haline stratification and increases in the baroclinic Rossby radius, making it harder for the shelf slope current to turn the tight corner into the Norwegian Trench. While we do find that the NWSPPE slope current tends to bypass the Norwegian Trench, we do not see the large-scale circulation configuration change in this region. We note that both Tinker et al. (2016) and Holt et al. (2018) found this while downscaling relatively coarse ocean models – Tinker et al. (2016) downscaled the 1.25° HadCM3 ocean and Holt et al. (2018) downscaled the 1° ORCA1 configuration of the NEMO model. Holt et al. (2018) also downscaled the $\mathrm{1}/\mathrm{4}$° ORCA25 configuration of NEMO and found this circulation change to be weaker and that the configuration did not flip to a new state. Our NWSPPE also downscaled ORCA025 (the ocean component of HadGEM3-GC3.05), so perhaps this northern North Sea configuration change, in part, depends on a coarser driving ocean model.

The wider North Atlantic shows a substantial freshening in the HadGEM3-GC3.05 PPE, which is consistent with the changes we see in the NWSPPE (Fig. S23). There are a number of mechanisms that could be contributing to this large-scale change, including changes to the hydrological cycle (Rhein et al., 2013), changes in stratification (Zika et al., 2018), Greenland Ice Sheet loss (Ruan et al., 2019) and AMOC slowdown (Zhu and Liu, 2020). The Atlantic freshening is consistent with the changes we see on the NWS. We have also shown that there is a slight reduction in the North Sea inflow (Fig. S12), which will also lead to a reduction of the NWS salinity. However, we do not see the North Sea circulation configuration change seen by Holt et al. (2018), which effectively isolates the North Sea from the Atlantic, with the North Sea becoming more estuarine. So while our slight reduction in the modelled NNSI will play a role in the NWS freshening, the main driver is likely to be the freshening of the North Atlantic.

As our NWSPPE is consistent with the PDCtrl, the two datasets can be used together. By using the projected climate evolution from NWSPPE, with the estimated unforced variability of the PDCtrl, it is possible, for example, to estimate the time of emergence of the climate signal from that climate variability. Most GCM climate projections are started from a convenient point within a climate spin-up simulation, when the modelled climate is stable (or its drift is acceptable), and the time counter is started from the year 1850 to reflect the start of the historical forcing period. Climate projections do not use data assimilation (or a similar process) to constrain the present-day climate to observed state of reality as used in seasonal or decadal forecasting systems (e.g. MacLachlan et al., 2014; Hermanson et al., 2022). This means that the present day of the climate simulation does not match reality, in terms of the phase of the various climate modes and variability. This does not matter in the distant future (late century) when we can be confident that the climate change signal dominates over the possible range of present-day variability, but this potentially limits the utility of climate projections in the near future due to the potentially confounding influence of internal variability (e.g. IPCC, 2021b). By estimating the time of the climate signal emergence (ToE), we can provide a lower bound of how soon you can start using uninitialised climate projections.

We have applied the method of Lyu et al. (2014) (described in Sect. 2.10) to give an estimate of the time of emergence of the SST and SSS climate signals from the present-day variability (Fig. S30), and we show the ensemble median and 16th–84th percentile range, masking out points where climate emerges in less than 84 % of the ensemble (10 of the 12 members). We find that the SST climate signal emerges from variability across the shelf with a relatively low (16th–84th) range across the ensemble, whereas the SSS signal only emerges from the climate variability for 42 % of the NWS (with emergence in greater than 10 ensemble members). Furthermore, the SSS climate signal emerges later than SST and with greater ensemble spread. To compare the ToE of NWS SST and SSS, we take the NWS means of the ToE median (and range) for the grid box ToE. We find that the SST typically emerges in 2034 (with an 8-year range), while SSS emerges in 2046 (with a 33-year range). This can guide where and when to use the NWSPPE for near-future temperature and salinity projections.

Table 9 Timings of exceeding global mean warming levels of 2 and 4 °C from the UKCP18 HadGEM3-GC3.05 PPE. The global warming levels are derived from the model simulated global annual mean anomaly relative to 1981–2000 baseline plus the observed warming from 1850–1900 mean to 1981–2000 mean based on HadCRUT4 (Morice et al., 2012). Timings are based on a centred 25-year running mean. While all simulations pass 2 °C of global mean warming, some simulations do not reach global mean warming levels of 4 °C by the end of the century. Adapted from Gohar et al. (2018).

UK policy-makers are interested in warming levels, such as the projected conditions under a 2 or 4 °C warmer world. These are a useful policy tool to describe the climatic conditions when the global mean temperatures rise to either 2 °C or 4 °C above pre-industrial levels. Warming levels have been used in the recent IPCC Fourth Assessment Report (IPCC, 2021b) and in the UK's third Climate Change Risk Assessment (CCRA; Betts et al., 2021) and have been calculated for the UKCP18 HadGEM3-GC3.05 PPE. Gohar et al. (2018) use a time-shifting methodology (Herger et al., 2015; Schleussner et al., 2016) to derive a pair of UKCP18-based 2 and 4 °C products. They use the HadGEM3-GC3.05 PPE (among other simulations) to find the years when the global mean temperature crosses the 2 and 4 °C thresholds relative to the pre-industrial conditions (considered to be 1851–1900). As we use the same PPE, the same years can be used as the basis of 2 and 4 °C projections of the NWS; however, as our simulations do not include the pre-industrial period, some care must be taken when trying to apply and interpret this methodology to the NWSPPE. Their Table 2 shows the timing of a centred 20-mean passing 2 and 44 °C of global mean temperature rise since the pre-industrial period for the anonymised models in their study (adapted in our Table 9). For example, their model 1 (our r001i1p00000) passes the 2 °C threshold in 2030 and the 4 °C threshold in 2063. While we can calculate the NWS conditions for r001i1p00000 in a 20-year period centred on 2030 (2021–2040) and can iterate through our ensemble, we cannot show how much the NWS has changed since the pre-industrial period, as our simulations start in 1980. Therefore, any differences between the NWS and the global mean that have occurred since the pre-industrial period are not captured. This is a common problem for regional studies. One approach is to say how different the NWS would be from the early century, under 2 and 4 °C warming levels. For this, for each ensemble member, we can use the time slices in Table 9 and subtract the 2000–2019 period and then average across the ensemble. Another approach is to add the global warming between 1850–1899 and 2000–2019 as calculated from the HadCRUT5 dataset (Morice et al., 2012), which equals 0.987 °C. As the NWSPPE does not begin in the pre-industrial period, the exact timings of the 2 and 4 °C warming levels will still remain uncertain, especially as the HadGEM3-GC3.05 PPE was not designed to span all sources of uncertainty. The exact methodology used to choose the period of the 2 and 4 °C warming levels may impact the averaging periods, but the methodology proposed here should give a first-order estimate.

These climate projections are a one-way forced downscaling of a GCM (HadGEM3-GC3.05) using climatological Baltic exchange and rivers. This approach is based on two main assumptions: firstly, that the driving GCM can simulate the large-scale climate and how it changes, and secondly, the downscaling approach allows for improved resolution and additional physical processes important to the shelf seas, without becoming inconsistent with the climate of the GCM.

HadGEM3-GC3.05 has been extensively evaluated at the global and regional scales (Williams et al., 2018) and for the parameters relevant for the NWS (Tinker et al., 2020). Dynamically downscaling GCMs for shelf sea regions is a well-established approach (Gröger et al., 2013; Mathis and Pohlmann, 2014; Tinker et al., 2016; Holt et al., 2010; Olbert et al., 2012). We follow the methodology of Tinker et al. (2020) using NEMO COx on AMM7, a very well-established model and domain. NEMO COx is used operationally at the Met Office for 6 d forecasts (O'Dea et al., 2017; Tonani et al., 2019) and in their reanalysis (Renshaw et al., 2019). Furthermore, it has been used as a research model on a number of timescales (including seasonal predictions, Tinker and Hermanson, 2021; centennial climate projections, Hermans et al., 2020b; and for the present-day unforced climate variability, Tinker et al., 2020). The choice of domain is important for downscaling – too big and the interior of the model will diverge from the climate of the parent model; too small and the boundary conditions will dominate the interior of the model. Furthermore, the geography of the domain is important; in the case of the NWS, shelf edge processes are complex and should be either excluded (e.g. Olbert et al., 2012) or be far enough from the LBCs to behave independently from them. As NEMO COx running on the AMM7 domain has been extensively used and evaluated, we are convinced that our approach is fit for our purposes, although we do note some limitations.

The exchange with the Baltic Sea is complex and helps control the dynamics of the Norwegian Trench, which is downstream of most of the NWS (including the UK). We use a climatology for this, which is a common approach (e.g. Tinker et al., 2015; Holt et al., 2010). The climatological Baltic may have an impact in the Skagerrak and along the coast of Norway, where the mean salinity (its variability) may be affected. This, with the use of climatological rivers, is an important limitation for understanding the NWS water cycle and salinity response as the climatological rivers and Baltic LBCs do not see any hydrological changes simulated in the parent GCM. This approach could be improved upon. A consistent set of transient Baltic Sea climate projections could give a trend to the climatology or even a set of transient boundary conditions. It may be possible to develop a Baltic Sea box model to simulate the exchange (e.g. Stigebrandt, 1987). The best approach would be to run both NWS and Baltic projections together. Care must be taken with the mean sea level of the Atlantic and Baltic LBCs, which can be inconsistent when taken from different models. We also use a river climatology, which is also a common approach, and is used operationally for the Met Office NWS forecasts (Tonani et al., 2019). HadGEM3-GC3.05 has online river routing, so using this river output directly may be a better option but would need assessment. The estimated Baltic outflow (668 ± 32 km³ yr⁻¹ for the period 2002–2021; Boulahia et al., 2022) has relatively low interannual variability and appears relatively small compared with the variability of the North Sea riverine inflow (NOSCCA, 2016) – this may suggest that improving the representation of the riverine climate response is a greater priority.

We use the HadGEM3-GC3.05 atmosphere directly. This has a spatial resolution of about 70 km, which is relatively low. Future projections may benefit from a higher-resolution atmosphere component; however, there may be little benefit to increasing the atmospheric resolution until you get to convection-permitting resolutions (∼ 2 km) (e.g. Vautard et al., 2013; Kotlarski et al., 2014). This could either come from using a higher-resolution GCM, output from a regional atmospheric climate model to downscale the atmosphere or by moving to a regional coupled climate model. Furthermore, other studies bias-correct the GCM model output before use (Mathis et al., 2013). This could be explored for uncoupled future studies.

We use COx on the AMM7 domain, which has a 7 km horizontal resolution, with 50 terrain-following vertical levels. AMM7 is an eddy-permitting model, but as the internal Rossby radius is ∼ 4 km on the NWS (Holt and Proctor, 2008), a model resolution greater than this would allow for the simulation of a much richer eddy field and improved representation of shelf exchange processes. AMM15 is a 1.5 km resolution model covering a similar region to AMM7. AMM15 will replace the use of AMM7 for many purposes. AMM15 is eddy resolving and shows an improved representation of the mean state across the domain compared with AMM7 (Graham et al., 2018). As well as having a higher resolution, it requires a smaller time step, so it is much more expensive to run computationally than AMM7. Currently, both AMM7 and AMM15 are run operationally at the Met Office to provide operational 6 d forecasts: AMM15 is coupled to a wave model, while AMM7 is coupled to the expensive biogeochemistry model, as AMM15 is currently too expensive to run with biogeochemistry. AMM15 may be too expensive to use as the basis for ensemble climate projections. However, it may be useful to run companion runs to show how these higher-resolution processes may influence the projections.

Many users are interested in estuaries or the coastal zone. These are currently beyond the scope of either AMM7 or AMM15. CO6 has a minimum depth of 10 m (in PDCtrl) to ensure the model does not dry out and crash. This is not necessary in CO9 (NWSPPE), where a wetting and drying module allows for more realistic coastal bathymetry. We suggest that locations within three grid boxes of the land are used with care and that the impact of the model's omission of coastal processes is considered on a case-by-case basis.

Our methodology uses an uncoupled, one-way forcing approach. While the HadGEM3-GC3.05 is coupled and so implicitly includes coupled processes in the COx model forcings, COx is uncoupled and so cannot feed back to the atmosphere or wider ocean. Coupled atmosphere–ocean processes can be very important in the tropics, and there is increasing attention to their importance in the NWS. For example, forecasting sea fog (Fallmann et al., 2019), ocean–wave coupling Lagrangian trajectories (Bruciaferri et al., 2021) and heat waves (Petch et al., 2020) can be improved by using a regional coupled model. The Met Office is developing a regional ocean–atmosphere-wave coupled model (UKC4) using the NEMO COx AMM15 model and the 2.2 km variable-resolution MOGREPS UK (Hagelin et al., 2017) with the RAL2 scientific options (Bush et al., 2023; Berthou et al., 2024). This will eventually couple with a biogeochemistry model and river model to allow for forecasting of multi-hazard compound events such as coastal flooding and erosion. This model is very complex and expensive to run, with the atmosphere component being much more expensive to run than the (already expensive) AMM15 model. While exploratory climate simulations of future time slices and case studies are possible, this is even less likely than AMM15 to be used as the basis of an ensemble of transient climate projections. However, its use with an ensemble of projections can help understand what processes are poorly represented and how the projections may overpredict or underpredict certain phenomena when compared with (things like) UKCP18.

Our NWSPPE gives a quantification of one source of climate uncertainty. By varying parameters within the atmosphere, land and aerosol components of HadGEM3-GC3.05, the range of the responses of the NWS climate can be assessed. This is an important source, but there are many other sources of uncertainty in climate projections. Other important sources include emission uncertainty, model structure uncertainty (of the global model, the shelf sea model and even the regional atmospheric climate model), initial condition uncertainty and methodological uncertainty (coupling and driving uncertainty). To date, the Minerva projections (Tinker et al., 2016) have been the only systematic and comprehensive assessment of climate uncertainty within the NWS (Tinker and Howes, 2020). While our NWSPPE updates the Minerva projections, they do not give any further insights into these other sources. A climate uncertainty budget for the NWS, which quantified and compared these sources, would be a useful contribution to the research field.

Our NWSPPE updates the Minerva projections by making use of more advanced models but does not extend their use. The next generation of NWS climate projections may couple to other components, such as a wave and/or a biogeochemistry models – both of which are coupled to COx for the Met Office operational synoptic forecasts. This will represent the (next) third-generation of NWS climate projections.

Our approach updates the Minerva NWS climate projections, which were based on those of Holt et al. (2010). We feel that these projections are needed for the UK's upcoming fourth CCRA, but that the next set of UK climate projections should include the third-generation NWS climate projections. These could address many of the limitations outlined above, by utilising improved resolution (AMM15), improved coupling methodology (UKC4), additional model systems (biogeochemistry, waves/surge) and different sources of uncertainty. These may require a range of methodologies, perhaps employing a cascade of model resolutions: AMM7 transient ensembles to address different sources of uncertainty, a transient simulation of the standard ensemble member with AMM15, and case studies and times slices with UKC4. While coupling to a wave model does not increase the computational expense much (∼ 10 %), adding biogeochemistry does (3–6 times), so biogeochemistry may be limited to the lower-resolution AMM7 (or with mixed resolution) simulations (as it is currently in operational forecasts). Adding biogeochemistry would allow a much wider range of end users to use the projections, while adding waves and surge would increase consistency and improve confidence. The next set of NWS climate projections may require a more community led approach to support their greater complexity. Close proactive engagement and collaboration with end users during the planning and development stage will enable uptake and the generation of user-relevant impacts to inform evidence-based coastal decision-making (Weeks et al., 2023).

7 Conclusions

Our key findings and conclusions are the following:

• The NWS annual mean SST rises by 3.11 °C (± 0.98 °C) for 2079–2098 relative to 2000–2019 and with a greater warming in the summer and autumn.

• The duration of the seasonal stratification increases by almost 1 month, initialising 4 d earlier and breaking down ∼ 19 d later.

• There is a region of limited warming in the deep water south of Iceland, associated with the AMOC slowdown.

• There is substantial freshening across the NWS and domain, with an NWS freshening of −1.01 (± 2σ = 0.93) for 2079–2098 relative to 2000–2019. The rate of freshening increases through the 21st century.

• There is a general reduction in strength of most of the NWS residual circulation over the 21st century and little change in its configuration, apart from the south and west of Ireland where there is a slight change in the configuration and a strengthening of the circulation; 80 % of the area to the south and west of Ireland shows an increase in the residual current magnitude, whereas 60 % of the remainder of the NWS shows a decrease.

• Using the NWSPPE with the PDCtrl allows the climate signal to be considered in the context of unforced natural variability. We give an example of how these can be used together to estimate the time of emergence of the climate signal from the natural variability. For example, the ToE for the annual mean SST averaged over the NWS is 2034, suggesting that the entire ensemble has moved out of the natural variability compared with 2000–2019.

• The ToE gives confidence to the use of our projections in the mid-century, and we discuss how they could be used to simulate changes consistent with 2 and 4 °C warming levels.

• These projections are for the physical environment, and so give no information about possible changes to the ecosystem (e.g. productivity, biomass, chlorophyll), biogeochemistry (e.g. nutrient levels, oxygen levels) or the wave conditions (e.g. significant wave height, wave period). In addition to including these other parameters, future marine climate projections may have greater spatial resolution, employ more sophisticated coupling strategies and assess additional sources of climate projection uncertainty.

Appendix A: Dataset description

The data underlying these NWSPPE projections, with the PDCtrl data, have been released on the CEDA website (https://catalogue.ceda.ac.uk/uuid/832677618370457f9e0a85da021c1312, last access: 20 May 2024). Here we will describe these files, the data structure and the available variables. For the NWSPPE, we provide monthly means of two-dimensional variables for every ensemble member, month, and year between 1990 and 2098. We also provide climatological means and standard deviations for each ensemble member for an early- and late-century period (2000–2019 and 2079–2098). We use these climatologies to provide ensemble statistics, which are the basis of many of our figures. We also provide regional mean time series.

NEMO is discretised onto an Arakawa C grid and so separates variables onto a T, U and V grid. Most scalar variables are on the T grid (the centre of the grid boxes). The U and V grids are offset between the T points (at the edges of the grid boxes). Therefore, variables on the T grid have a different location (in terms of longitude and latitude) compared with their equivalent U and V grids. We use this distinction and have three sets of files, with most variables in the T grid files, and only the U and V components of the barotropic current in the U and V grid files. We expect the T grid files to be sufficient for most purposes and applications. We also calculate the barotropic current magnitude on the T grid (reported in the T grid files) – this gives users a first look at the circulation.

A1 Directory structure

Here we describe the directory structure. At the highest level, the NWSPPE data (Tinker, 2023b) are separated from their resulting ensemble statistics (Tinker, 2023a) and the PDCtrl data (Tinker, 2023c). Under the NWSPPE data, there are sub-directories for each ensemble member, each of which have sub-directories for the annual files (annual), the regional mean files (regmean) and the climatologies (clim). PDCtrl has sub-directories for the annual files (annual) the regional mean files (regmean):

NWSClim

1) NWSPPE

   a) r001i1p00000

      i) annual

      ii) regmean

      iii) clim

(repeated for all 12 NWSPPE ensemble members)

   b) …

2) EnsStats

3) PDCtrl

   a) annual

   b) regmean

A2 File names

Here we give example file names to help explain our naming convention.

NWSPPE/r001i1p00000/annual

• NWSClim_NWSPPE_r001i1p00000_1990_gridT.nc

• NWSClim_NWSPPE_r001i1p00000_1990_gridU.nc

• NWSClim_NWSPPE_r001i1p00000_1990_gridV.nc

These are for the T, U and V grid files, respectively, where the ensemble number and year change.

NWSClim/NWSPPE/r001i1p00000/regmean

• NWSClim_NWSPPE_r001i1p00000_1990-2098_regmean.nc

These are for the regional mean files. The regional mean files give monthly means for the whole period.

NWSClim/NWSPPE/r001i1p00000/clim

• NWSClim_NWSPPE_r001i1p00000_clim_xxx_YYYY-YYYY_gridT_mean.nc

• NWSClim_NWSPPE_r001i1p00000_clim_xxx_YYYY-YYYY_gridT_stddev.nc

These are for the climatological mean and standard deviation, respectively, on the T grid – again the U and V grid variables are in separate files. The xxx denotes the monthly (m01-m12 for January to December); seasonal (djf (December to February), mam (March to May), jja (July to August), son (September to November) for winter, spring, summer or autumn, respectively); or annual means (ann), and the YYYY-YYYY denotes the climatological period and may be 2000–2019 or 2079–2098.

We note that the variable names are unchanged (SST in one file is the 20-year mean, and in the other it is the 20-year standard deviation). The long_name variable attribute is also unchanged, but the method is captured in the cell_method.

NWSClim/PDCtrl/annual

• NWSClim_PDCtrl_1990_gridT.nc

• NWSClim_PDCtrl_1990_gridU.nc

• NWSClim_PDCtrl_1990_gridV.nc

The PDCtrl monthly mean files are analogous to the NWSPPE files described above.

NWSClim/PDCtrl/regmean

• NWSClim_PDCtrl_2050-2250_regmean.nc

The NWSPPE regional mean files are analogous to the NWSPPE files described above.

NWSClim/EnsStats

• NWSClim_NWSPPE_EnsStats_clim_xxx_YYYYYY
  YY_gridT_stats.nc

For the ensemble statistics (on the T grid, U and V grid variables are in separate files), where xxx again denotes monthly, seasonal or annual means. While the ensemble statistics are produced for the same climatological period (i.e. YYYYYYYY can be 2000–2019 or 2079–2098), they are also produced for the difference between the periods, so YYYYYYYY can also be 2079-2098minus2000-2019.

Here the file contains four statistics per variable, and so the variable names are modified to give the statistic name – SST_ensmean, SST_ensvar, SST_intvar and SST_ensstd. The long name and cell methods also reflect the statistic.

A3 Variables

The variables that we have released include sea surface (SS-) and near-bottom (NB-) temperature (T) and salinity (S) and their difference (DF-) (SST, NBT, DFT, SSS, NBS, DFS), potential energy anomaly (PEA), mixed layer depth (MLD) and sea surface height (SSH), the U and V components of the barotropic (depth mean) velocity (DMU and DMV), and barotropic (depth mean) speed (DMUV). These are all on the T grid, apart from DMU and DMV on the U and V grids, respectively.

The regional mean time series include RegAveSST, RegAveNBT, RegAveDFT, RegAveSSS, RegAveNBS, RegAveDFS, RegAveSSH and RegAvePEA (the regional means of SST, NBT, DFT, SSS, NBS, DFS, SSH and PEA, respectively) – all on the T grid variables, apart from DMUV and MLD.

Additional variables and frequencies may be available to collaborators; contact the corresponding author for details.

A4 File structure

Here we describe the monthly mean fields, climatologies and ensemble statistics available at the CEDA data centre.

The released dataset covers the period between 1990 and December 2098 and is in CF-compliant NetCDF file format (CF-convention 1.8).

A4.1 Annual (two-dimensional fields)

Monthly mean data is available for all months from January 1990 to December 2098, for all 12 ensemble members, and the present-day control simulation from 2050–2250 (each year still representing the year 2000). For each year, there are three files for the model T, U and V grid, with 12 monthly means within each. All variables are on the T grid, apart from the U and V components of the barotropic velocity (on the U and V grid, respectively). Only 2D data fields are being released, with 3D variables such as temperature and salinity being reduced to surface (SS-) and near-bed fields (NB-), as well as the difference between them (DF-).

Details of the files are given in Table A1 in terms of the NetCDF variable names and their dimensions. The dimensions of the files are “time”, “lon”, “lat” and “bnds” (simplified to t, y, x and 2, respectively, in Table A1). Time is given in seconds since 1950-01-01 00:00:00 UTC, using a 360 d calendar. The time variable (time) is supplemented by the time bounds (time_bounds), giving the start and end of the averaging period (hence having two dimensions). The 2D data fields are references with the longitude (lon) and latitude (lat) variables in each of the files – these are slightly offset between the T, U and V grids, although this is unlikely to cause any difficulty in most applications. Time is given in as double-precision (“f8”) floating-point numbers, whereas other variables are given as single precision (“f4”) floating-point numbers.

A4.2 Climatologies

The monthly mean files are used to produce 20-year early-century (2000–2019) and late-century (2079–2098) climatologies, giving the mean and standard deviation in separate files. These are produced for annual means (ann), seasonal means (djf, mam, jja and son) and monthly means (m01–m12 for January to December). The NetCDF files have the same format as outlined in Table A1. The time_bounds variable sets the start and end of the climatology. For climatologies with an even number of years, the average time can be misleading, especially for monthly or seasonal climatologies where, for example, a winter climatology time_counter will be for summer. However, if you think of 2 sequential winters, the average time (in seconds since 00:00 UTC 1 January 1950) will be the summer between them – this is an unfortunate, if correct, feature of climatologies with an even number of years.

A4.3 Ensemble statistics

These climatologies are used to produce the ensemble mean, ensemble variance, and standard deviation and interannual variance. As the NWSPPE is designed to capture the likely uncertainty (spread) associated with uncertain parameters rather than to give the most likely outcome, the ensemble mean is not the most likely value. However, the ensemble mean and the ensemble standard deviation are a useful way to summarise the (typically Gaussian) distributions. While the ensemble standard deviation is simply the square root of the ensemble variance, we think it will be a more useful statistic, so we include it to save end users from having to calculate it. We also give the interannual variability, following Tinker et al. (2016). These data are given for the near-present-day (2000–2019) and future (2079–2098) periods. For the difference between the periods (2079-2098minus2000-2019), the statistics are simply the difference between the future and present-day statistics. This is useful to show how the ensemble has changed (i.e. how the ensemble mean and ensemble standard deviation have changed); however, this does not give information on the uncertainty associated with a projected change. We therefore provide two additional statistics in the ensemble statistics difference files: the projected ensemble mean (projensmean) and the projected ensemble standard deviations. Here, we remove the present-day climatological mean from the future climatological mean to give a resulting anomaly ensemble – we then calculate the ensemble mean and standard deviation. This is useful when we want to say how the SST has increased, with an estimate of the uncertainty on the projection; for example, when we say the southern North Sea winter SST increased by 3.55 °C (± 1.20 °C), we are taking winter southern North Sea regional mean of SST_projensmean (3.55 °C) and (2×) SST_projensmean (1.20 °C).

We consider the late-century climate change by comparing the early century (2000–2019) to the late century (2079–2098) for monthly, seasonal and annual means. These time slices have been updated since UKCP09 and the Minerva projections, partly to make them closer to the present but also as the simulations started in 1980 rather than 1952.

Following Tinker et al. (2016), we use 2 times the standard deviation (2σ) to describe the ensemble spread (average ± x, where x = 2σ). We calculate 2σ for each grid box for a given month (and season and year) across each of the 20-year time periods and plot these in Figs. 9–13. The typical interannual variability, however, is often greater than the ensemble time mean spread. Tinker et al. (2016) developed a more sophisticated analysis of the ensemble variability by decomposing the total variance (${\mathit{\sigma }}_{\text{tot}}^{\mathrm{2}}$) into an interannual variability component (${\mathit{\sigma }}_{\text{int}}^{\mathrm{2}}$), which is the mean of the 20-year interannual variances for each ensemble member, and an ensemble spread component (${\mathit{\sigma }}_{\text{ens}}^{\mathrm{2}}$), which is the variance of the 20-year means for the ensemble:

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\mathit{\sigma }}_{\text{int}}^{\mathrm{2}}={\displaystyle \frac{{\sum }_{e,y}\left((x_{e,y}{)}^{\mathrm{2}}\right)}{n_e{n}_y}}-{\displaystyle \frac{{\sum }_e\left({\left(\frac{{\sum }_y{x}_{e,y}}{n_y}\right)}^{\mathrm{2}}\right)}{n_e}},\text{(A2)}& {\displaystyle }{\mathit{\sigma }}_{\text{ens}}^{\mathrm{2}}={\displaystyle \frac{{\sum }_e\left({\left(\frac{{\sum }_y{x}_{e,y}}{n_y}\right)}^{\mathrm{2}}\right)}{n_e}}-{\left({\displaystyle \frac{{\sum }_{e,y}{x}_{e,y}}{n_e{n}_y}}\right)}^{\mathrm{2}},\text{(A3)}& {\displaystyle }\begin{array}{rl}{\mathit{\sigma }}_{\text{tot}}^{\mathrm{2}}=& {\mathit{\sigma }}_{\text{int}}^{\mathrm{2}}+{\mathit{\sigma }}_{\text{ens}}^{\mathrm{2}}\\ =& {\displaystyle \frac{{\sum }_{e,y}\left((x_{e,y}{)}^{\mathrm{2}}\right)}{n_e{n}_y}}-{\displaystyle \frac{{\sum }_e\left({\left(\frac{{\sum }_y{x}_{e,y}}{n_y}\right)}^{\mathrm{2}}\right)}{n_e}}\\ & +{\displaystyle \frac{{\sum }_e\left({\left(\frac{{\sum }_y{x}_{e,y}}{n_y}\right)}^{\mathrm{2}}\right)}{n_e}}-{\left({\displaystyle \frac{{\sum }_{e,y}{x}_{e,y}}{n_e{n}_y}}\right)}^{\mathrm{2}}\\ =& {\displaystyle \frac{{\sum }_{e,y}\left((x_{e,y}{)}^{\mathrm{2}}\right)}{n_e{n}_y}}-{\left({\displaystyle \frac{{\sum }_{e,y}{x}_{e,y}}{n_e{n}_y}}\right)}^{\mathrm{2}},\end{array}\end{array}$$

where this decomposition uses the sum of squares formula and where x_e,y denotes a variable (such as SST) for a given year y and ensemble member e; n_y is the number of years in the sample (20), and n_e is the number of ensemble members (12).

For each variable in the monthly means, there are now three ensemble statistics, so rather than produce three files we have added a suffix to the variable name.

A5 Regional means

Maps of climatological means are good at showing spatial patterns and changes but do not show how the system evolves temporally. Regional mean time series show how the system evolves, show how the ensemble spread changes and also give an estimate of the NWS response to model parameter uncertainty. They also allow for a first look at the data to help inform a more detailed analysis. We follow the methodology of Tinker et al. (2019) by averaging the model fields over a relevant region mask every time step. We use the region mask of Wakelin et al. (2012) as it divides the NWS into regions that make geographic and oceanographic sense. We also include a “shelf” region, by combining several NWS regions (excluding the Atlantic, Norwegian Trench, Skagerrak/Kattegat and Armorican shelf regions) – when we refer to a NWS mean or averaged over the NWS, we refer to this shelf region. NEMO calculates these regional means every time step and then averages them into the monthly mean time series (hourly means are also available), which are released for SST, NBT, DFT, SSS, NBS, DFS, PEA and SSH. We also include the Wakelin et al. (2012) region mask in the file (mask). The time, time_bounds, lat and lon (for the mask) variables are also in other files. The variables are prefixed by RegAve (e.g. RegAveSST, RegAveSSS, RegAvePEA) and have two dimensions – time and region. Two one-dimensional variables reg_id and cnt specific to region mean files are included. reg_id is the region ID and is a value from 0 to 13 – this aligns with the region names (given in the global attribute region_names). cnt is a count of the number of grid boxes within each region (ranging from 792–25 523).

One difference in the regional means from NEMO 3.6 and NEMO 4.0.4 (used in PDCtrl and NWSPPE, respectively) is that NEMO3.6 simply averaged the values of all the grid boxes in a region, while NEMO 4.0.4 takes an area-weighted average. For a direct comparison between the two datasets, it is fairly easy to reprocess the regional means from the annual files, with or without area weighting. Example code is provided in the GitHub Python package.

A6 Sea surface height (SSH)

Care must be taken when using the projected SSH – this is not to be used directly as a sea level projection. Sea level is affected by a wide range of processes, and so sea level projections are built up from several different models (Fox-Kemper et al., 2021). Furthermore, given the different resolutions between HadGEM3-GC3.05 and COx, the different sources of lateral boundaries for the Atlantic and Baltic, and the absences of tides in the GCM, the SSH in COx is likely to diverge from HadGEM3-GC3.05. However, COx improves the representation of coastal ocean processes, and so it does give additional information on sea level change (e.g. Hermans et al., 2020a, b; Tinker et al., 2020), but this is only part of the sea level puzzle (Palmer et al., 2018).

The SSH data provided here are obtained from downscaling the HadGEM3-GC3.05 simulation of ocean dynamic sea-level (DSL). DSL is defined as the local height of the sea surface above the geoid and is corrected to remove time-dependent variations of the sea surface due to atmospheric loading referred to as the inverse barometer effect. In the driving model (HadGEM3-GC3.05), DSL is defined to have zero global ocean area mean; therefore, in the downscaled simulations the SSH data are providing the local anomaly in geocentric sea level that results from local ocean density (steric) and circulation (dynamic) effects. While COx does simulate tides and the inverse barometer effect, its representation of SSH change includes only part of the drivers of sea-level change active in the real world, which we briefly outline below.

Sea level change arises from processes that alter the total volume of water in the ocean column or processes that alter the mass of water by redistributing water masses between the oceans and the continental land-surface. Simulated changes to ocean density from climate models can be used to estimate the effect of changes in ocean volume. However, many important processes that alter ocean mass, notably the contributions from loss of land-ice mass are not included in climate models, and estimates for these process must be supplied from offline models. Additionally, changes to the distribution of water mass in the oceans and over the land surface will produce spatially varying patterns of sea-level change associated with the earth's gravitational field, axial rotation and deformation of the solid-earth surface. The effect of the gravitation, rotation and deformation process, referred to collectively as GRD, must also be estimated using offline models. Finally, local sea-level projections for a coastal location may refer to changes relative to a local solid surface rather (relative sea-level change) than a fixed geoid (geocentric sea-level change) and use additional models to account for local vertical movement of the land surface (e.g. subsidence). These disparate sources are pulled together into the sea level puzzle (Palmer et al., 2018).

A7 Missing data

Over the 110 years, 12 months and 12 ensemble members, there were 23 files that were not archived correctly. Most of these were recreated by averaging the daily mean files giving near-identical files (with the difference at the noise level). Three files (Table A2) could not be created this way, and so they were replaced by averaging the monthly mean for the same month from the previous and following years. The approach should maintain the seasonal cycle and any trends in the data; however, it may impact some statistics, such as estimates of interannual variability. Full details of the replacement files are given in Tables A3 and A4.

A8 Calculation of the barotropic current speed on the T grid (DMUV)

The use of model currents can often be complicated, especially as the U and V grids are offset. Often maps of current speed are useful, as they give an illustration of circulation. We have therefore provided the barotropic current magnitude on the T grid. First, we interpolated the U and V components of the barotropic current on the T grid:

$$\begin{array}{}\text{(A4)}& \begin{array}{r}\mathrm{DMU}\mathit{\_}{T}_{j,i}={\displaystyle \frac{{\mathrm{DMU}}_{j,(i-i)}+{\mathrm{DMU}}_{j,i}}{\mathrm{2}}},\\ \mathrm{DMV}\mathit{\_}{T}_{j,i}={\displaystyle \frac{{\mathrm{DMV}}_{j,i}+{\mathrm{DMV}}_{(j-\mathrm{1}),i}}{\mathrm{2}}},\end{array}\end{array}$$

where DMU and DMV are the U and V components of the barotropic current on their native U and V grids. Once they have been converted to the T grid, they become DMU_T and DMU_T, respectively. We then calculated their magnitude (DMUV):

$$\begin{array}{}\text{(A5)}& {\mathrm{DMUV}}_{j,i}=\sqrt{\mathrm{DMU}\mathit{\_}{T}_{j,i}^{\mathrm{2}}+{\mathrm{DMU}\mathrm{\_}\mathrm{T}}_{j,i}^{\mathrm{2}}}.\end{array}$$

A9 Code

This post-processing was undertaken in Python and is available from GitHub:

https://github.com/hadjt/NWS_simulations_postproc (last access: 20 May 2024)

A10 Digital object identifier (DOI)

Each of the top-level datasets has been issued with a DOI.

Physical Marine Climate Projections for the North West European Shelf Seas: NWSPPE. 20 July 2023.

https://doi.org/10.5285/edf66239c70c426e9e9f19da1ac8ba87 (Tinker, 2023b).

Physical Marine Climate Projections for the North West European Shelf Seas: PDCtrl. 20 July 2023.

https://doi.org/10.5285/66e39885a60e4b6386752b1a295f268a (Tinker, 2023c).

Physical Marine Climate Projections for the North West European Shelf Seas: EnsStats. 20 July 2023.

https://doi.org/10.5285/bd375134bd8c4990a1e9eb6d199cc723 (Tinker, 2023a).

A11 Appendix A tables

Table A1 Details of NetCDF files.

Table A2 Files missing from the NWSPPE that we recreated by averaging adjacent files.

Table A3 Monthly mean two-dimensional files that were not archived and how they were replaced: “Daily Means” shows that the daily mean files were archived and have been averaged into monthly means.

Table A4 Regional mean monthly mean files not archived.

Appendix B: Tides in 360 d calendar

The dynamics of the NWS and other shelf seas are dominated by tides. They provide the mixing energy to seasonally mix the water column in winter and define whether a region is mixed or stratified in the summer (Simpson and Bowers, 1981). They drive important aspects of the residual and instantaneous circulation (Tinker and Hermanson, 2021). They interact with the sediment and coastline, leading to geomorphological processes (Pingree and Griffiths, 1979).

Tides are controlled by aspects of the orbits of the earth–moon–sun systems (Pugh, 1987). The tidal frequencies are based on five astronomical frequencies, including the tropical year, the sidereal day and the sidereal month (Pugh, 1987). Different combinations of these frequencies give the frequencies of the tidal constituents, such as M₂ and S₂, that dominate in most regions (Pugh, 1987). The addition of these terms, and the beating between them, leads to important behaviours, including the spring–neap cycle and the equinoctial tides (Pugh, 1987).

Climate models use a variety of calendars. In addition to the Gregorian (365 plus leap years) calendar, other common calendars include the 365 d with no leap year calendar and the 360 d calendar. The 360 d calendar is made up of twelve 30 d months which make post-processing and calculation of monthly and seasonal means much simpler. However, the 360 d calendar is not directly compatible with the astronomical constants of the real world. This is not very important from the climate modellers point of view, as most current climate models do not have tides – from a global climate perspective, realistic, dynamic tides are not particularly important, and their role can generally be performed by a mixing parameterisation (Tinker et al., 2022). It becomes a problem when climate models are dynamically downscaled for shelf sea regions (e.g. Holt et al., 2010; Tinker et al., 2016). As global climate models advance, tides will become increasingly integrated; then this may also become a problem with GCMs.

When running a shelf sea model with a 360 d calendar, the default approach (used within NEMO) is to reset the tides at the beginning of every month, so the equinoctial tides always occur at the correct time of the year (hereinafter “reset”; see Table B3). However, this introduces a jump in the tides as days are repeated or skipped (for non-30 d months), which could introduce noise into the system and renders tidal analysis (based on least-square multiple linear regression) useless. Another approach is to convert the 360 d model forcings into 365 d and to make use of the standard Gregorian calendar. This is typically achieved by repeating days or interpolating (through time) to stretch the forcings (hereinafter “Greg”; see Table B3). Both have drawbacks, including affecting the speed of synoptic systems. Another approach is to run the tides continuously from the start of the run (on their Gregorian calendar) and accept that they will become out of sync with the model running on the 360 d calendar, gaining ∼ 5 d every year (hereinafter “drift”; see Table B3). This means that the equinoctial tides drift through the year. A final approach is to adjust the astronomical frequencies and so the tidal constituent frequencies to match a 360 d year, which correctly fixes the equinoctial tides to a point in the year (hereinafter “compress”; see Table B3). Each of these approaches have strengths and weaknesses, which may influence the simulations. Here we describe the three tidal implementations (reset, drift, and compress) and their effects on the tidal characteristics of the NWS. We do not assess the modification of the atmospheric and oceanic forcings.

We first give a brief overview of relevant aspects of tidal theory: how the astronomical constants may be adjusted for the 360 d year and the impact of this on the tidal constituent periods and the important beating frequencies. We then assess the impact of this tide in a model simulation and compare it to the other methodologies.

For the NWSPPE and PDCtrl, we take a conservative approach and use the “drift” method, following Holt et al. (2010) and Tinker et al. (2016).

B1 Method

B1.1 Tidal theory

Here we give a brief overview of relevant aspects of tidal theory, mainly taken from Pugh (1987).

Harmonic expansion of the tide breaks it into many sinusoidal constituents with astronomically defined frequencies and amplitudes and phases that vary with space. For a given location, the effect of the tide on the sea level can be defined:

$$\begin{array}{}\text{(B1)}& \mathit{\eta }=\sum a_n\sin ({\mathit{\sigma }}_{n}t+{\mathit{\varphi }}_n)+{\mathit{\eta }}_{\mathrm{0}},\end{array}$$

where η is the sea level (and η₀ is the mean sea level), and for constituent n, a_n is the amplitude, σ_n is the frequency ($\mathrm{2}\mathit{\pi }/T$ where T is the tidal period) and φ_n is the phase. The a_n and φ_n can be found with least-square multiple linear regression (Pugh, 1987).

There are a number of important astronomical frequencies that define the frequencies of the tidal constituents: the mean solar day (${\mathit{\omega }}_{\mathrm{0}}=\mathrm{2}\mathit{\pi }/\mathrm{1}\mathrm{d}=\mathrm{2}\mathit{\pi }/\mathrm{24}\mathrm{h}$), the mean lunar day (${\mathit{\omega }}_{\mathrm{1}}=\mathrm{2}\mathit{\pi }/\mathrm{1.0351}\mathrm{d}=\mathrm{2}\mathit{\pi }/\mathrm{24.8424}\mathrm{h}$), the sidereal month (${\mathit{\omega }}_{\mathrm{2}}=\mathrm{2}\mathit{\pi }/\mathrm{27.3217}\mathrm{d}=\mathrm{2}\mathit{\pi }/\mathrm{655.7208}$), the tropical year (${\mathit{\omega }}_{\mathrm{3}}=\mathrm{2}\mathit{\pi }/\mathrm{365.2422}\mathrm{d}=\mathrm{2}\mathit{\pi }/\mathrm{8765.8128}\mathrm{h}$), moon's perigee (${\mathit{\omega }}_{\mathrm{4}}=\mathrm{2}\mathit{\pi }/\mathrm{8.85}$ years), the regression of the moon's nodes (${\mathit{\omega }}_{\mathrm{5}}=\mathrm{2}\mathit{\pi }/\mathrm{18.61}$ years) and the perihelion (${\mathit{\omega }}_{\mathrm{6}}=\mathrm{2}\mathit{\pi }/\mathrm{20}\mathrm{942}$ years).

The tidal constituent frequencies (ω_n) can be calculated:

$$\begin{array}{}\text{(B2)}& {\mathit{\omega }}_n=i_{\mathrm{1}}{\mathit{\omega }}_{\mathrm{1}}+i_{\mathrm{2}}{\mathit{\omega }}_{\mathrm{2}}+i_{\mathrm{3}}{\mathit{\omega }}_{\mathrm{3}}+i_{\mathrm{4}}{\mathit{\omega }}_{\mathrm{4}}+i_{\mathrm{5}}{\mathit{\omega }}_{\mathrm{5}}+i_{\mathrm{6}}\mathit{\pi }{\mathit{\omega }}_{\mathrm{6}},\end{array}$$

where i₁, i₂, i₃… are small integers (|i| < 5). The use of both the solar and lunar days is not necessary as ${\mathit{\omega }}_{\mathrm{0}}={\mathit{\omega }}_{\mathrm{1}}+{\mathit{\omega }}_{\mathrm{2}}-{\mathit{\omega }}_{\mathrm{3}}$. The species of the constituent is defined by i₁ (i₁=2 for semi-diurnal, 1 for diurnal, or 0 for long period) (Pugh, 1987). Doodson notation concisely describes the tidal constituents: i₂−i₆ are incremented by 5 (so there are no negative integers); i₁−i₆ are combined into a single number. For example, M₂ has a Doodson number of 255.555 (i₁=2, ω_n=2ω₁), S₂ is 273.555 (i₁=2, i₂=2, $i_{\mathrm{3}}=-\mathrm{2}$, ${\mathit{\omega }}_n=\mathrm{2}{\mathit{\omega }}_{\mathrm{1}}+\mathrm{2}{\mathit{\omega }}_{\mathrm{2}}-\mathrm{2}{\mathit{\omega }}_{\mathrm{3}}$) and K₂ is 275.555 (i₁=2, i₂=2, ${\mathit{\omega }}_n=\mathrm{2}{\mathit{\omega }}_{\mathrm{1}}+\mathrm{2}{\mathit{\omega }}_{\mathrm{2}}$),

When two tidal constituents of slightly different frequencies are added together, they beat as they come in and out of phase. The M₂ and S₂ constituents beat to give the spring–neap cycle, and the S₂ and K₂ constituents beat to give the equinoctial tides. The beat envelope function is

$$\begin{array}{}\text{(B3)}& \pm \left[a_{\mathrm{1}}+a_{\mathrm{2}}\left[\cos \left(\right({\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{2}})t+({\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}\left)\right)\right]\right],\end{array}$$

where the constituents are ordered by amplitude size (i.e. a₁ > a₂). This gives the beat frequency (σ₁−σ₂) and phase (φ₁−φ₂), so the M₂–S₂ spring–neap cycle has a frequency of $\mathrm{2}\mathit{\pi }/\mathrm{12.42}-\mathrm{2}\mathit{\pi }/\mathrm{12}=\mathrm{354.85}\mathrm{h}=\mathrm{14.7}\mathrm{d}$, and the S₂–K₂ has a beat frequency of 182.6221 d.

You can also calculate the beat frequencies from the Doodson numbers. For example, the difference between M₂ and S₂ is 255.555, and 273.555 implies a difference in frequency (${\mathit{\omega }}_n=\mathrm{2}{\mathit{\omega }}_{\mathrm{2}}-\mathrm{2}{\mathit{\omega }}_{\mathrm{3}}$), while S₂ and K₂ is 273.555, and 275.555 is ${\mathit{\omega }}_n=-\mathrm{2}{\mathit{\omega }}_{\mathrm{3}}$.

B1.2 Compress tidal scheme

We introduce a tidal scheme where one of the astronomical frequencies is modified to the 360 d year. The frequency of the tropical year (ω₃) is modified to $\mathrm{2}\mathit{\pi }/\mathrm{360}\mathrm{d}=\mathrm{2}\mathit{\pi }/\mathrm{8640}\mathrm{h}$. This has no impact on the S₁ or S₂ tidal frequencies but reduces the M₂ period by 1.84 s. The impact on the beat frequencies is larger, with the spring–neap cycle increasing from 14 d 18 h 22 min to 14 d 18 h 47 min and the equinoctial tidal frequencies decreasing from 182.621 to 180 d, exactly half a 360 d year.

We have implemented this scheme in the shelf sea model NEMO Coastal Ocean model 9 (CO9, NEMO version 4.0.4), and we will assess the impact of this simulation of the tides. Note that we have not used this scheme in the NWSPPE.

B1.3 Models and experimental design

We use the NEMO version 4.04, Coastal Ocean model version 9 (Wise et al., 2024), run on the 7 km NWS domain (Atlantic Margin Model 7 km – AMM7) in a set of 40-year simulations with different approaches to modelling the tide (see Table B5). We use these simulations to investigate the impact of the different tidal schemes on the resulting emergent tidal behaviour. We assess these runs for their amphidromic systems and distributions of amplitude and phase of the different tidal constituents. We also pay particular attention to the difference in amplitude and phase between key constituent pairs, as this controls the behaviour of their tidal beats, which is a particular weakness of some of the tidal schemes We do not assess the impact on other oceanic properties (such as temperature), which is likely to be more subtle, and require longer simulations to investigate.

We assess the second 20-year period used in this assessment, as this ties in with the inline tidal harmonic analysis within NEMO (within DIA/diaharm_fast.F90, within https://github.com/hadjt/NEMO_4.0.4_CO9_shelf_climate/tree/master/src/OCE/, last access: 20 May 2024).

B1.4 NEMO model code

Most of the code added to NEMO is within SBC/tide_mod.F90 (https://github.com/hadjt/NEMO_4.0.4_CO9_shelf_climate/blob/master/src/OCE/SBC/tide_mod.F90, last access: 20 May 2024). This allows the user to choose between the three tidal schemes using the nam_tides360 namelist and the ln_tide_drift and ln_tide_compress logical keywords.

B1.5 Analysis techniques: least-square tidal harmonic analysis

We use a tidal harmonic analysis based on a least-square method to calculate the phase and amplitudes of the tidal constituents as described by Pugh (1987). This method considers the tidal signal as

$$\begin{array}{}\text{(B4)}& T\left(t\right)=Z_{\mathrm{0}}+\sum _N{H}_{\mathrm{n}}{f}_n\cos [{\mathit{\sigma }}_{n}t-g_n+(V_n-u_n\left)\right],\end{array}$$

where the Z₀ is the mean sea level, H_n and f_n are the amplitude and the nodal factor for the nth tidal constituent, t is time σ_n is the tidal frequency g_n is the tidal phase, and V_n and u_n are the phase angle at time zero and the nodal angle, respectively. t and the terms f_n, σ_n, V_n, and u_n are known a priori, as is the surface elevation T(t). Equation (B4) is solved for Z₀, H_n, and g_n by using a least-square estimation.

This is done online by NEMO (see https://github.com/hadjt/NEMO_4.0.4_CO9_shelf_climate/blob/master/src/OCE/DIA/diaharm_fast.F90, last access: 20 May 2024) using consecutive 20-year periods of the model simulation. We use the second 20-year period to avoid any spin-up issues.

One issue of the standard NEMO tide configuration (“Reset”) is that the time counter jumps at the end of 360 d calendar months for months that should have 31 d (or February). NEMO's online harmonic analysis code takes account of this, as both the tide, and the time counter jump. We have written Python code that replicate the NEMO online harmonic analysis code. We can then see how these jumps affect the tidal analysis if not taken into account – i.e. if someone simply analysed the hourly sea surface height model output. This is used in Fig. B2c (cf. Fig. B2b), which shows very large RMSE, as well as (not shown) a near-zero M₂ amplitude and larger S₂ amplitude.

B2 Results

B2.1 Harmonic analysis

We fit a least-square tidal fit to the data (Fig. B1). We use the RMSE of the residuals to quantify how well the data is fit (Fig. B2). As the tides see discontinuous time for Reset, with days repeated and missing, we carefully recreate the time array for its tidal analysis; otherwise, the most diurnal constituents (including M₂) appear to be S₂. We include the RMSE when this is not done (see Fig. B2c), which shows an ∼ 10-fold increase in RMSE. Generally, all models have similar RMSE, which are small compared with the tidal range, and so we can use the computed phase and amplitudes for the rest of the study.

The locations of the amphidromes are insensitive to the simulation (Fig. B3). The difference of the amplitudes between the simulation is typically less than 5 %. The phases of the Compress experiment have the largest difference.

The phase of individual diurnal and semi-diurnal constituents is not thought to be of any particular importance as the peak will move through the day (apart from S₁, S₂ – hence spring high water is always at the same time). However, the phase difference between constituents give the phase of the tidal beats, which can be important, and can occur at the same time of the year.

B2.2 Important beat frequency and amplitude

To identify the important beat frequency pairs, we use Eq. (B3) to find the beat frequency of each pair of constituents and their amplitude (a₂ in Eq. B3, i.e. the smaller of the two amplitudes) for the shelf mean (Fig. B4). We focus on beat pairs with periods greater than 100 d. Figure B4 shows that the S₂–K₂ pair are the dominant biannual beat pair, with K₁–P₁ playing a secondary role. The annual beat frequencies are much weaker, with the S₂–T₂ pair dominating and K₁–S₁ and P₁–S₁ having a lessor (and similar) role. K₂–T₂ is tri-annual, while N₂–NU₂ and 2N₂–MU₂ occur every 202 d 10 h 20 min and 202 d 13 h 35 min, respectively.

Figure B4 and Table B6 show that the proposed “Compress” tidal system correctly adjusts the beat frequencies to the 360 d year; however, the phase of the beat is also important – it is what makes the equinoctial tides occur during the equinoxes. We can use Eq. (B3) to find the beat phase (the difference between the individual tidal constituent phases) and calculate the timing of the peak of the beats – we use this to show how well the different tidal systems reproduce the phase of the beating waves.

Figure B5 shows that the Gregorian and the Reset systems have the same timing of the maxima. The Drift can be seen to drift throughout the period, progressing ∼ 5 d yr⁻¹. Compress is comparable to Greg and Reset, showing that it is capable of reproducing the period and timing of the tidal beats. We can look at the maps of the spatial patterns of the timing of the beat pair maximum (Fig. B6). For the most important pair (S₂–K₂), the patterns of the timing agree visually and the timing is within ∼ 10 d, which is considered reasonable.

B3 Conclusions

We have compared the ability of different 360 d tidal schemes to simulate the emergent behaviour of the NWS tides from simulation with the standard Gregorian calendar. We base our analysis on the results of a least-square tidal harmonic analysis. This has to be modified for two tidal schemes (using the modified tidal periods for Compress and using a synthetic time for Reset).

Reset is the standard tidal scheme within NEMO but leads to jumps within the tide at the end of the month. This changes the average M₂ period, and so the errors associated with the (uncorrected) harmonic analysis are much higher (Fig. B2c).

With the corrected the harmonic analysis, all 360 d tidal schemes simulate the cotidal charts of the control (Greg) well (Fig. B1) with similar (pattern and magnitude of) RMS errors (Fig. B2). This is the most important quality, and each scheme performs equally well.

The beat phases give the timings of the equinoctial tide. Both Reset and Compress simulate the correct timings, whereas the beat period does not match the 360 d calendar. So even if the phase is correct, the timings drift through the calendar each year, typically gaining 5 d every year.

For both these criteria, the Compress tidal scheme performs well, and is not out-competed by any other scheme. However, it is unpublished, and some users may be impacted by some of its features (M₂ tide shortened by 2 s, spring–neap cycle increases by 25 min). We therefore use the Drift tidal scheme in our projections, as it has been used in earlier projections (e.g. Holt et al., 2010; Tinker et al., 2016) and is a more conservative choice.

B4 Appendix B tables

Table B1 The 360 d tide methodologies and their advantages and disadvantages.

Table B2 Period of updated tidal constituent

Table B3 Name and description of simulations.

Table B4 Table B4 periods of beat pairs

B5 Appendix B figures

Figure B1 Cotidal charts for the M₂ (a–d), S₂ (e–h), K₂ (i–l) and S₁ (m–p) tidal constituents under the Gregorian (control, a, e, i and m), Reset (b, f, j and n), Drift (c, g, k and o) and Compress (d, h, l and p) tidal schemes. The colour maps represent the amplitude and have constant colour limits for a given constituent, while the lines are lines of equal phase (every 45°).

Figure B2 RMS of least-square tidal fit for the simulations: (a) Gregorian; (b) Reset; (c) Reset (uncorrected); (d) Drift; (e) Compress. For Reset, we construct an equivalent time (resetting to the Gregorian calendar at the beginning of each month) and use this for the least square fit. We do not correct this in Reset (uncorrected), and the RMS errors are much greater (c). Note that the colour scales vary between panels.

Figure B3 Cotidal charts for M₂, S₂ and K₂ (a–d, e–h and i–l, respectively) for Gregorian (a, e, i) and anomalies (relative to Gregorian) for Reset (b, f, j), Drift (c, g, k) and Compress (d,h, l). The colours give the amplitude (or amplitude minus Gregorian amplitude anomaly) in mm; the contours give the co-phase lines (every 45°): grey for Gregorian and white for Reset, Drift, and Compress, respectively.

Figure B4 Beat frequency and amplitude for the standard Gregorian tidal constituent frequencies as averaged over the NWS (a), and the modified 360 d frequencies under the Compress tidal scheme (b). The yellow and red lines show the monthly and bimonthly frequencies, and the blue and black lines show the biennial and annual frequencies ($\mathrm{2}\mathit{\pi }/(\mathrm{365.24}\times \mathrm{24})$ and $\mathrm{2}\mathit{\pi }/(\mathrm{360}\times \mathrm{24})$, respectively). The grey box shows the beat pairs with periods less than 100 d and amplitudes less than 0.2 m.

Figure B5 The timing (day of the year) of the maximum of the beating tidal constituent pairs (with periods greater than 100 d) for the four different tidal systems for an example location. Each plot is a beating tidal pair. The black line represents the day of the year (in the x direction) and the decimal year (in the y direction). The red dot is the local maximum of the beating pair. Each tidal system (Greg, Reset, Reset (uncorrected), Drift, and Compress) is stacked on top of one another.

Figure B6 Timing of largest beat pair maximum (in day of the year), for pairs that occurs at approximately the same time of each year. Note that S₂-K₂ are the beating pair with the greatest amplitude.