1 Introduction

Knowledge about sea surface height variability in the form of probability density functions (PDFs) is an important ingredient for understanding and projecting sea level change under global warming conditions (Church et al. 2013; Stammer et al. 2019). Information about complete sea level PDFs is especially useful to stakeholders and decision makers as it contains information about potential future high-end sea level changes, information that is embedded in the upper tails of distribution functions and resides outside the likely range usually provided by the Intergovernmental Panel on Climate Change (IPCC) (Stammer et al. 2019). High-end sea level changes pose great risk for flooding and highest costs in adaptation efforts (Le Cozannet et al. 2017; Hinkel et al. 2019); yet, respective information is incomplete, due to a hitherto missing data basis required to infer full sea level PDFs that account for all processes involved. Consequently, there exist little consensus on changes in regional probabilistic sea level, and on potential consequences for coastal risk management (Rasmussen et al. 2018; Oppenheimer et al. 2019). In fact, not much is known about the impact of climate change on sea level statistics as only limited statistical sea level information exists even for today albeit respective regional observational and modeling studies (Johansson et al. 2001; Dangendorf et al. 2017; Lang and Mikolajewicz 2019, 2020).

In the absence of any detailed PDF tail-information, high-end sea level estimates are so far based either on probabilistic approaches (Kopp et al. 2014; Jevrejeva et al. 2019) or on expert solicitations (Bamber et al. 2019). Faced with shortcomings intrinsic to both, Stammer et al. (2019) introduced a new conceptual approach for estimating future high-end sea level rise beyond the likely range as function of time scale and emission scenarios considered by practitioners according to their risk-averseness. Building on this approach, van de Wal et al. (2022) recently provided first estimates of high-end numbers as function of time scale. Despite of this progress, strong demands remain for an actual estimate of PDFs of sea level variability and for information about their changes in a warming world. In particular, understanding the extent to which future extreme sea level rise will occur still requires an improved knowledge of regional probabilistic sea level projections that provide the full information on possible sea level variability. We recall in this context that at any moment high-end sea level (extreme) events depend on slow changes of the climatological state due to climate change superimposed by natural variability which will amplify the extremes. This paper is about the potential changes of the natural variability contribution to extreme sea level events.

Estimating probablistic sea level variability information typically involves computing PDFs of sea surface height (SSH) and its changes (changes in distribution of likely range and high-end tails called percentiles) from either observations (Sura and Gille 2010) or Coupled Model Intercomparison Project (CMIP)-type climate models. Contemporary probabilistic sea level variability information, estimated previously from satellite altimetry (Sura and Gille 2010), shows substantial geographic variations in sea level variability statistics in particular in eddy-rich regions. Concerning future regional sea level change patterns, Slangen et al. (2014) analyzed sea-level changes based on CMIP5 scenarios amended by off-line information on contributions from the cryosphere and the solid Earth. Respective analyses are also available along coastlines providing estimates of sea level distribution functions and coastal relative sea-level (RSL) change statistics (Carson et al. 2016) that could support of coastal adaptation considerations (Haasnoot et al. 2019).

Recently improved approaches toward more reliable regional sea level statistics were made possible through the initiation of large ensemble simulations (Deser et al. 2012; Becker et al. 2023) allowing to estimate full sea level distribution functions in the model world as well as to detect significant changes in SSH statistics in space and time (Widlansky et al. 2015; Frankcombe et al. 2018). Today, several large ensemble simulations exist generated using single coupled models by applying small changes to the initial conditions of each member to establish the full ensemble spread (Maher et al. 2019; Deser et al. 2020). These ensembles were run for historical conditions and under various climate change radiative forcing scenarios. Building on previous variability research involving large ensemble climate model simulations, here we aim to use the results to further improve our knowledge about sterodynamic sea level statistics as a function of geographic location, thereby also linking changes in statistics to adjustments of underlying dynamical processes.

In detail, we use the Max Planck Institute Grand Ensemble (MPI-GE) simulations (Maher et al. 2019) to estimate monthly changes in sea level probability density functions (PDFs) and their regional variations from the 100 ensemble members, as well as related statistical properties, such as variance, skewness and kurtosis, the latter two being indicative of the Gaussian nature or deviations thereof. Deviations from Gaussian nature measured by skewness and kurtosis can help gain insights of the mechanisms behind driving ocean circulation (e.g., Sura and Gille 2010) and how these mechanisms might change under climate forcing conditions. To this end, we investigate changes in projected regional sea level statistics (variability) by the end of the 21st century simulated under the two emission scenarios RCP4.5 and RCP8.5, such as shifts in skewness and kurtosis, as well as the probability distribution statistics at selected locations. But we note that sea level variability originating from any process not included in the model will not be accounted for, such as ocean eddies, storm surges, changes in ocean tides, or sea level contributions originating from the solid earth in terms of sea floor uplift or gravity changes, nor contributions arising from terrestrial ice mass loss or terrestrial water storage changes. Hence, we do not claim to cover all aspects of the evolution of extreme sea-levels under global warming but focus on the changing characteristics of sterodynamic sea level.

The remaining paper is structured as follows: in Sect. 2 we describe the methodology, including the MPI-GE model and the statistical analysis method of the model output. Changes of sterodynamic sea level statistics under climate change conditions are discussed in Sect. 3. Mechanisms behind changes of sea level variability in the tropical Pacific, North Atlantic and the Southern Ocean are described in Sect. 4. Concluding remarks follow in Sect. 5.

2 Method

2.1 The MPI grand ensemble

The Max Planck Institute Grand Ensemble (MPI-GE; Maher et al. (2019) is based on a variant of the MPI Earth System Model (MPI-ESM), Version 1.1 (Giorgetta et al. 2013), used during the fifth phase of the Climate Model Intercomparison Project (CMIP5) but in detail is closer to the MPI-ESM Version 1.2 used during CMIP6 (Mauritsen et al. 2019). The MPI-GE ocean component is implemented at a horizontal resolution of 1.5° and 40 vertical layers with grid poles located over Greenland and Antarctica. Given this resolution, the model is not capable of representing ocean eddies. The atmosphere model ECHAM6 has a spectral horizontal resolution of T63 (equivalent to ~ 1.8° spatial resolution) and 47 vertical layers. The MPI-ESM 1.1 has been evaluated previously against observational records for the atmosphere (Giorgetta et al. 2013), the ocean (Jungclaus et al. 2013) and sea ice (Notz et al. 2013).

As part of the MPI-GE effort, experiments were carried out using historical forcing conditions as well as RCP4.5 and RCP8.5 climate-change forcing, each experiment encompassing 100 ensemble members. Historical and RCP simulations follow the experimental design guidelines of CMIP5 (Taylor et al. 2012). This MPI-GE model output has been analyzed previously to assess responses in Northern Hemisphere stratosphere after tropical volcanic eruptions (Bittner et al. 2016), the influence of internal variability on global radiation budget (Hedemann et al. 2017) and changes in Artic sea-ice in a 1.5 °C warmer world (Niederdrenk and Notz 2018). Here we use the output to study sea level statistics in its largest extent possible hitherto.

To compute sea level statistics from the GE output, we use the model’s sea surface height (SSH) field in its sterodynamic form defined as the sum of the global mean steric and the dynamical sea surface height on the native model grid (Gregory et al. 2019). Changes in statistical properties of sterodynamic sea level were assessed as a function of climate forcing by examining differences in statistical parameters as they emerge between the historical period (1986–2005) and the end of the two RCP 4.5 and 8.5 projections (2081–2100). Since initial states of ensemble members are sampled from the control simulation (Table 1 of Maher et al. 2019), model drift introduces substantial offsets between the members leading to artificially enhanced ensemble spread. To account for this model drift, we estimated and removed the mean off-set and linear trend from the pre-industrial control simulation from each ensemble member start year (taken from the list of start years from Table 1 of Maher et al. 2019). In this way, the influence of the drift of pre-industrial control simulations on the sea level statistics was removed from each ensemble member (of historical, RCP4.5 and RCP8.5 simulations) prior to computing statistical parameters.

2.2 Statistical analysis

Taking advantage of the large ensemble, throughout this paper, statistics are computed on the ensemble dimension, i.e., at each month over all ensemble members. Computations include (1) changes in ensemble mean and standard deviation (SD) (2) 95th and 99th percentiles relative to the ensemble mean, (3) timeseries analysis, (4) ensemble spread as well as global spatial maps for SSH statistics over 20 years of between the historic period (1986–2005) and the RCP projections (2081–2100). To infer percentiles and higher statistical moments (skewness and kurtosis), probability distribution functions (PDFs) are computed in the form of histograms on the combined time and ensemble dimension to capture the total variability. For this purpose, occurrence frequencies were computed based on counts of regional SSH anomalies. The bin-size is calculated based on the maximum and minimum in the dataset of all three distributions (historical, RCP4.5 and RCP8.5) to create an array of equally spaced bins. This results in the same number of bins with their bin width changing according to how the whole PDF width changes and hence consistently representing the same percentage of PDF per bin for all selected regions and periods. For each PDF, the three central statistical moments were evaluated: (i) variance (ii) skewness (S) and (iii) excess kurtosis (K’). Dimensions of variance are in cm^2; S and K’ are dimensionless. Positive values of S imply that the PDF is skewed to the right, i.e., toward high monthly sea level, the opposite holds for negative skewed values. Following in Bulmer (1979), we define S with amplitudes between ± 0.5 as “approximately symmetric”, between ± 0.5 and ± 1.0 as “moderately skewed” and between > + 1.0 or < −  1.0 as “highly skewed”. The excess kurtosis (K’) is defined as K-3, with 3 corresponding to the kurtosis for a Gaussian distribution. A Gaussian (symmetrical) distribution would have zero S and K’ of 0. For K > 3, the PDF is more peaked than a Gaussian distribution, whereas for K < 3, it is flattened.

During the following analyses, the statistical significance of fields of the mean sterodynamic sea level, and its higher moments SD, S and K’ and their changes will be determined at the 95% confidence level using a bootstrap method with 5000 iterations of randomly selected monthly anomalies having the same size as each group of the mean, SD, S, K’ by month. For each field and at any given grid point, the monthly mean, SD, S, K’ differences between the historical and RCP runs are identified as significant at the 95% confidence level if changes lay outside the 2.5–97.5 percentile range of the bootstrap distribution. To evaluate if the distributions for the two periods are significantly different from each other, we used Student’s t-test (to test if the sample means are different) and an f-test (to test whether the sample variances are from different populations/distributions) thereby following Zwiers and von Storch (1995) for the calculation of the equivalent sample sizes. From the PDFs, we also show the shifts in the high-end tails in the 95th percentiles.

3 Regional SSH statistics and its changes

To set the stage for the later analysis of statistical SSH parameters and their changes under global warming conditions, shown in Fig. 1 are the ensemble mean fields of the open ocean sterodynamic SSH averaged over the 20-year period (2081 to 2100) for the RCP 4.5 and 8.5 runs, respectively, minus the 20-year mean from the historical run (1986 to 2005). In comparison to the CMIP5 multi-model ensemble mean SSH change for RCP4.5 (see Church et al. (2013); their Fig. 13.13a), MPI-GE shows a comparable regional pattern suggesting a decreased subtropical gyre circulation in the North Atlantic (Gulf Stream region) and an increase in the North Pacific (Kuroshio region) with similar patterns seen in eddy evolution for both regions (Beech et al. 2022). However, sea level anomalies are more pronounced in our single-model (MPI-GE) ensemble mean which agrees to what can be observed for sterodynamic change resulting from individual CMIP5 models (see Church et al. (2013); their Fig. 13.18).

Fig. 1
figure 1

Difference between a RCP4.5 and b RCP8.5 (2081–2100), respectively, and the historic run (1986–2005). The global mean values are ~ 0.21 m in the top panel and ~ 0.35 m in the lower panel for RCP 4.5 and 8.5, respectively. The averaged ratio of the two panels amounts to a factor of 1.5. White grid points mark anomalies that are statistically insignificant at the 95% confidence level.

Full size image

Ensemble mean maps of SSH standard deviation (SD) are shown in the top row of Fig. 2, as they emerge under (a) historical, (b) RCP4.5 and (c) RCP8.5 conditions. Shown in the middle (Fig. 2d–f) and bottom rows (Fig. 2g, h) are fields of S, and K’, for each of the above three 20-year periods. For Gaussian statistics to hold, fields of S and K’ should be zero. Clearly this is approximately the case for most of the global ocean. However, it does not hold in the tropics, especially in the Pacific, where S < − 1 in western tropical area and S > 1 in the east. With respect to K’, a clear deviation from its Gaussian value of zero can be found in the western tropical Pacific, exceeding 3.5 there as it is the case for PDFs broader than bell-shaped normal distributions. As pointed out by Sura and Gille (2010), for stochastically forced systems, non-Gaussian statistics could be characteristic for the non-linearity of the system.

Fig. 2
figure 2

a–c Present and future sterodynamic sea level SD (in cm) computed from mean monthly sea surface height anomalies over the historical (1986–2006), and from the RCP4.5 and RCP8.5 runs, respectively, representing the period (2081–2100). d–f S and g–i K’ for each grid box of mean monthly sea surface height anomalies for the historical period (1986–2006), RCP4.5 (2081–2100), RCP8.5 (2081–2100). Color bar increments in the bottom row differ for positive and negative values. Fields of S and K’ do not have units. Six regions of interest highlighted as black squares are selected and analyzed as PDFs (see Fig. 5)

Full size image

To provide a sense of the ensemble spread ability to reproduce sea surface height statistics observed by altimetry in the time domain, we compare the model results with the statistics available from Copernicus Climate Change Service (C3S) altimeter data for the period 1993–2013 for 20 years of sea level statistics (Fig. 3). The C3S data set consist of gridded monthly mean sea level anomalies which are computed with respect to a 20-year mean reference period (1993–2012) at ¼ degree resolution. The standard deviations seen from Fig. 3a, b confirm that the MPI-GE model roughly provides the same spatial pattern as the observed SSH variability; however, the model underestimates the amplitudes of open ocean SSH variability, especially along the western boundary current regions due to the lag of ocean eddies in the coarse-resolution model. Overall, the MPI-GE variance amplitudes therefore remain at a 25% level of what is being observed as our model is not capable of simulating eddies. Hence, all statistical considerations presented in this paper are thus confined to large-scale variability patterns and conclusions drawn about changes of extreme sea level events are likely to be an underestimate of what will happen in reality.

Regarding the model performance in estimating patterns for S, Fig. 3c, d suggest an overall agreement in the tropics with respect to the asymmetric pattern across the tropical Pacific. This is also in agreement with other studies (e.g. Sura and Gille 2010) which link this with ENSO oscillations, leading to negative tails in the west and positive tails of the distributions in the east due to strong El Niño—La Niña asymmetry in amplitude. Beech et al. (2022) suggested that regions of high kurtosis coincide with regions of high kinetic energy seen in eddy-rich regions. The inability of the MPI-GE in reproducing the observed patterns of K’ is therefore expected.

The changes of the statistical parameters under RCP4.5 and RCP8.5 conditions relative to the historical run are shown in Fig. 4. In both projections, the variability increases relative to the historical run over large parts of the Atlantic, the Arctic, the Pacific and the Indian Ocean (Fig. 4a–c). However, it decreases in the western North Atlantic and the Pacific and Atlantic sectors of the Southern Ocean (Fig. 4a–c). In the difference between RCP8.5 and the historical run the variability also decreases over the western tropical Pacific Ocean. In contrast, variability increased markedly in the RCP8.5 run southwest of Australia; a fact that also shows up between the two RCP runs, albeit smaller in amplitude. Compared with the historical period, RCP8.5 shows the largest reduction in both S and K’ for the western tropical Pacific region, which can be linked with changes in future ENSO activities (see discussion Sect. 4).

Fig. 3
figure 3

Mean Monthly sea-level statistics between the C3S sea level anomalies and ensemble statistics for the MPI-GE historical SSH for the period 1993–2013 for SD (a, b), S (c, d) and K’ (e, f). The altimeter data [taken from the Copernicus Climate Change Service (C3S)] consists of inter-calibration process of Topex/Poseidon between 1993 and 2022, Jason-1 between 2002 and 2008, Jason-2 between 2008 and 2016

Full size image
Fig. 4
figure 4

Difference in statistics between (left) RCP 4.5, and (right) RCP8.5 minus historical period. a–c Standard deviation (SD) (in cm) of mean monthly sea surface height anomalies over the historical (1986–2006), and from the RCP4.5 and RCP8.5 runs, respectively representing the period (2081–2100). d–f S and g–i K’ for each grid box of mean monthly sea surface height anomalies for the historical period (1986–2006), RCP4.5 (2081–2100), RCP8.5 (2081–2100). Grid points are masked in white if the anomalies are not statistically significant at the 95% confidence level

Full size image

Fields of S and K’ shown in Fig. 2 led us to conclude that generally the model sea level variability for this 20 years’ time scale appears to be close to Gaussian. This holds for the historical period (1986–2005) as well as the RCP4.5 and RCP8.5 (2081–2100) projections. However, a few regions stick out as non-Gaussian as highlighted by the regions shown in Fig. 2. To evaluate this non-Gaussian nature of the SSH variability in these regions in more detail, we show in Fig. 5 estimates of the underlying probability density function (PDF) computed as histograms of the SSH anomalies in several dynamically distinct regions marked as black box in Fig. 2. Respective statistics are summarized in Table 1.

Table 1 Changes in sea level PDF statistics related with Figure 5
Full size table

The figure confirms our earlier findings of variance changes resulting in different shapes of the PDFs (seen from shifts between the historical, RCP4.5 and RCP8.5 curves). It appears that the impact of global warming results in widened tails for the eastern North Atlantic and narrower PDF with contracted tails in the western North Atlantic of the RCP8.5 compared to the historical PDF (Fig. 5a, b). This widening in the PDF tails is especially true for the subtropical North Pacific (Fig. 5e), which appears to be linked to the presence of enhanced sea level variability. In contrast, seen for the western Pacific and the South Pacific regions (Fig. 5c, f), S and K’ in Table 1 do not deviate from Gaussian characteristics under global warming conditions. For the western Pacific, the historical period ENSO related SSH variability leads to PDFs highly skewed to the left but this behavior reverses in the RCP8.5 scenario (see S and K’ in Table 1), making the PDF more Gaussian as discussed above. The opposite pattern is seen for the eastern Pacific (Fig. 5d). The changes seen here for the PDFs agree with the Fig. 4 patterns and can be linked with similar physical mechanisms.

Fig. 5
figure 5

Regional SSH PDFs from six locations selected from Fig. 4d of a Northwest Atlantic, b Northeast Atlantic, c Western Pacific, d Eastern Pacific, e North subtropical Pacific and f South Pacific area. The vertical lines represent the 5th and 95th percentiles of the respective distribution

Full size image

Statistics associated with the PDFs shown in Fig. 6. The statistical significance of change in variance between periods was assessed using an f-test at the 95% level where f-test1 is between the RCP 4.5 and historical, and f-test 2 is between RCP 8.5 and historical.

Changes in open ocean SSH statistics under global warming conditions are typically assumed to be small and changes in high-end monthly sea level are therefore proposed to be linked just with changes in the mean SSH (Menéndez and Woodworth 2010). Under these assumptions, future changes in regional mean SSH have been used to predict the frequency of the occurrence of extreme events simply by shifting the relation between extreme sea level and return period derived from extreme value distributions upward by the expected mean sea level rise (Hunter 2010, 2012; Lowe and Gregory 2010; Lin et al. 2016; Garner et al. 2017). Using the MPI-GE results, we can now test the hypothesis of whether the statistics of sea level variability will change or not under global warming conditions. To this end we have defined extreme monthly changes by selecting two different percentiles (the 95th and the 99th ) and estimating the change of each percentile relative to the ensemble mean for both RCP runs. While the selection of both percentiles (95th and 99th ) may appear somewhat arbitrary, we believe that these upper ends of the PDF can prove useful in extracting changes in open ocean monthly high-end sterodynamic sea level changes linked to large scale mechanisms which can serve as vital information for coastal adaptions under climate change.

Shown in Fig. 6a, b is the respective differences relative to the mean for open ocean SSH PDFs of the historical run. Overall, the figure reproduces the variability pattern over the ocean indicating that the regional differences in (95th ) 99th percentiles is governed primarily by differences in SSH variability rather than other changes of the PDF characteristics. The figure suggests that under present-day conditions, natural variability enhances the likelihood of high-end sea level to occur in high-variability regions, especially along western boundary currents of the ocean.

To quantify future shifts in the high-end tails under climate warming conditions, Fig. 6b–d show the same diagnostics from the last 20-year period of the RCP4.5 and RCP8.5 runs relative to the historical SSH 95th and 99th percentile. All four panels reveal that high-end tails contract in a warmer climate over the Gulf Stream along the east coast of the U.S. and Canada, confirming findings from Fig. 4b of a reduced variability there. Tails also contract over parts of the Southern Ocean, e.g., the Ross Sea. In contrast, RCP8.5 changes suggest that most of the global ocean will see increased tail values in the future due to an increase of SSH variability under global warming conditions. This is most obvious for the Southern Ocean southwest of Australia, the eastern North Atlantic as well as the Kuroshio/Oyashio region.

Due to lack of statistical information, the assumption of unchanged statistics is often made and only the mean is adjusted to explain consequences of warming. With this assumption, in a graph that plots the return period vs. sea level extreme, one only needs to shift the curve along the axis of the amplitude of sea level extreme. Our results for the high-end changes seen as the percentiles-mean, show that a significant portion of the changes might be related to a change in variability and not just a shift in future mean changes. In this context, we recall that the underlying model is only coarse resolution and does not produce ocean eddies. Its variability therefore is substantially underestimated as compared to altimetric SSH observations (see also below). However, when using eddy-resolving model simulations in the CMIP6 framework, a recent study by Beech et al. (2022) provides a mixed evaluation of long-term evolution of ocean eddies in a warming world; with eddies increasing in some regions and decreasing in other regions (see Beech et al. (2022); their Fig. 2). Moreover, important processes such as storm-surges or tide are not included in the model. Hence, future changes in the extreme sterodynamic sea level tails might therefore be more severe in the real world than suggested here.

Fig. 6
figure 6

a Map of the SSH 99th percentile relative to the ensemble mean for the last 20 years of the historical run (1986–2005), b the same but for 95th percentiles—mean. Shown in c, d and e, f are similar fields for RCP4.5 and RCP8.5 representing the period (2081–2100) from which each of the fields shown in (a, b) were subtracted

Full size image

4 Sea level variability change mechanism

Figures 4 and 6 suggests that several geographic regions show pronounced changes in SSH variability under climate change conditions, and in its Gaussian characteristics. Notably, these are the tropical Pacific, the North Atlantic and the Southern Ocean, south-west of Australia. In the following we will analyze in more detail the variability changes in these regions.

4.1 Tropical Pacific

In the tropical Pacific, much of the open ocean regional sea level variability has been attributed to internal variability in the climate system associated with the El Niño Southern Oscillation (ENSO) (Stammer et al. 2013; Han et al. 2017). This link between SSH variability and ENSO was also identified previously by Sura and Gille (2010) and Widlansky et al. (2015, 2020). ENSO is known to dominate this region by altering surface winds, ocean currents, increased sea surface temperature (SST) variability, salinity and wind-driven shifts of the thermocline, which in turn can affect SSH (Widlansky et al. 2015, 2020). As shown in Maher et al. (2018, 2022), and in Ng et al. (2021), the MPI-GE run does show a clear ENSO contribution to variability in the tropical Pacific. We can therefore hypothesize here that any key changes in SSH statistics can hint at changes ENSO activity. But it is also worth noting that the magnitude of contributions from other climate modes, such as the wind driven Pacific Decadal Oscillation also affect the Pacific basin sea level (see Han et al. (2017).

As pulses of the easterly wind stress are among the important mechanisms for initiating ENSO events, we investigate changes in the statistics of the easterly wind stress to occur in Fig. 7. We find a consistent change in zonal winds stress statistics, which suggests that the S originating from easterly wind stress pulses reduces under RCP8.5 warming conditions. Similarly, the reduction in the historical positive K’ for both SSH and wind stress also points to a reduction of high-end monthly changes for both quantities. Our analysis suggests that a possible reduction of ENSO activity in the western tropical Pacific, especially under RCP8.5 forcing. However, it could be plausible that future ENSO intensity does not weaken altogether as influences from other mechanisms may influence parts of the Pacific region, e.g., for increasing central Pacific (CP) El Niño events as suggested by modeling studies (Trenberth and Stepaniak 2001; Ashok et al. 2007; Kug et al. 2009). Cai et al. (2018) diagnose a robust increase in the RCP8.5 El-Niño SST variability and suggested this increase is largely due to greenhouse gas-warming-induced intensification in the equatorial Pacific. However, large single model ensembles (Maher et al. 2019; Ng et al. 2021) and multi-model ensembles under different IPCC scenarios (Cai et al. 2014, 2018, 2022) suggests that the range in ENSO variability is large and differs widely between model simulations.

Fig. 7
figure 7

Differences in statistics of mean monthly surface easterly wind stress showing the difference between the RCPs minus historical period, as well as RCP8.5 - RCP4.5. a–c Standard deviation (SD) (in cm) d–f S and g–i K’ for each grid box. Grid points masked in white indicate anomalies that statistically are not significant

Full size image

4.2 North Atlantic

While searching for dynamical causes for the large-scale changes in variance in the North Atlantic, we compare in Fig. 8 the standard deviations of SSH with that of SST. Visually, both rows reveal a similar resemblance of changes in variance, with decreases in western North Atlantic, and increases in eastern North Atlantic, suggesting that both changes are linked to changes in circulation pattern associated with climate modes. The resemblance holds especially for the western North Atlantic, where the decrease in variability may be related to a shift in the North Atlantic current system eastward. The increase in SSH variability in eastern North Atlantic can be linked with future increases in the variability of surface ocean warming.

Fig. 8
figure 8

Changes in SSH standard deviations (cm) plotted as difference between the variability during RCP runs minus those obtained from the historical period, as well as the respective difference between RCP8.5-RCP4.5 between (a–c). d–f same but now as SST standard deviation patterns SST (in °C). Grid points are masked in white if the anomalies are not statistically significant at the 95% confidence level.

Full size image

Han et al. (2017) suggests that North Atlantic basin-scale SSH variability is driven by seasonal and interannual processes that dominate the North Atlantic Gyre circulation. Notably this would entail coupled atmospheric-ocean patterns, such as the North Atlantic Oscillation (NAO) (Hurrell and Deser 2009; Deser et al. 2017), evolution of ocean eddies (Beech et al. 2022) and changes between wind and heat fluxes driven variations in the Atlantic Meridional Overturning Circulation (AMOC) (Häkkinen 2001; Han et al. 2017, 2019). We also note that under global warming conditions, the MPI-GE projects a decrease in both the ensemble mean AMOC strength at 26°N (compared to 19 Sv for the historical run in 2000, the RCP8.5 declines to 14 Sv, while the RCP4.5 declines to 16 Sv by the end of 2100) and, similarly, its variability (ensemble spread) over time (see Maher et al. (2019); their Fig. 10). Since AMOC and the gyre circulation are dynamically coupled in the subpolar North Atlantic (Yeager et al. 2015), changes in both strength and variability characteristics of AMOC will translate to the strength and variability of the subpolar gyre.

Several studies have also highlighted that increased heat and freshwater fluxes in the subpolar North Atlantic weaken the AMOC and cause a reduction of AMOC strength and northward ocean heat transport (Jungclaus et al. 2014; Lobelle et al. 2020). The decreased sea levels seen in the western North Atlantic can be associated with a reduced transport of the Gulf Stream and the North Atlantic Current that leads to declining sea level gradients across these currents with a large scale impact on sea level (Häkkinen 2001; Brunnabend et al. 2014; Landerer et al. 2014; Hu and Bates 2018; Hand et al. 2020). A future declining AMOC and weakened westerly winds (Fig. 9a–c) can potentially be responsible for weak/slow circulation and poleward shift of the Gulf Stream.

First EOF pattern shows no significant change in the MPI-GE NAO strength under both scenarios (not shown). However, using the ensemble members from the CESM1 climate model outputs, for the RCP4.5 and RCP8.5, Hu and Bates (2018), show no change in the pattern, magnitude, or position of the NAO but relate the reduced intra-ensemble sea level variability to the decrease in AMOC variability. Moreover, future changes in sea level pressure and the North Atlantic jet Stream (westerly wind stress) patterns related with the NAO are predicted to change by shifting poleward and weakening (Deser et al. 2012, 2017; Woollings and Blackburn 2012; Hu and Bates 2018). This, in turn, can influence the gyre circulation and the structure of the AMOC and can influence the North Atlantic (increase-decrease) SSH variations. Also, Beech et al. (2022) predict that future eddy activities intensify poleward for the North Atlantic (shown in their Fig. 2c) and along most eddy rich areas like the Western Boundary Currents. However, this includes a mixed outlook, with both increases and decreases shown regionally.

4.3 Southern Ocean

The most pronounced changes in future SSH variability in the RCP8.5 scenario occurs in the Southern Ocean south-west (SW) of Australia [48–56°S, 78–123°E]. In this region, the SSH variability increases more rapidly since mid 1960s and undiminished until 2100 for RCP8.5 while the increase saturates in the mid-century for the RCP4.5 runs (Fig. 9). This variance increase is primarily related to changes of interannual variability and tapers off towards shorter and longer time scales, suggesting that interannual changes in the flow field become very prominent in this region characterized by strong topographic interaction with the Antarctic Circumpolar Current (ACC) flow field.

Previously, several studies using CMIP5 models and satellite altimeter have linked changes in the large meandering mesoscale band south of Australia [100−140°E] to shifts in the southern frontal position of the ACC transport as well as variability in the Southern Annular Mode (SAM) (Gille 1994, 2014; Landerer et al. 2014). Regarding specifically the region SW of Australia [100–140°E], several studies have linked changes in the large meandering mesoscale band to shifts in the southern ACC front (Gille 1994, 2014; Landerer et al. 2014) as well as the strengthening and southward shift in westerly wind stress on sea level change (Bouttes et al. 2012). The SSH variance increase in the MPI-GE under RCP8.5 forcing SW of Australia may be linked to increases in projected zonal eastward wind stress (see Figs. 7a–c, 9a–c and 10). Moreover, altimetry and reanalysis data also shows that SSH variability is linked to both the Southern Oscillation and the SAM (Armitage et al. 2018) and seem to be reflected in the southward shifts in southern ACC front and associated topography interactions (Sokolov and Rintoul 2009).

The large scale changes based on CMIP5 models suggest that changes of the Southern Ocean circulation are driven by a combination of factors such as wind stress variability, changes in horizontal circulation in form of the mean path and strength of the ACC, and associated interactions between topography and poleward eddy transports in the subpolar gyre area (Meijers 2014). This is shown by recent quantification of historical and projected Southern Hemisphere westerly wind changes based on analysis of CMIP5, CMIP6 and reanalysis data (ERA 5) by Goyal et al. (2021), suggests CMIP6 models have less equatorward bias in the westerlies as compared to CMIP5, and an associated reduction in projected poleward shift when compared to CMIP5 models (their Fig. 1a and b). Also, under the higher emission scenarios (e.g., SSPS-8.5) the CMIP6 ensemble mean shows a stronger jet strength compared with the CMIP5 ensemble mean (see Goyal et al. 2021; their Fig. 1b) and this poleward intensification projects onto a positive trend in the SAM (their Fig. S3). Other studies such as Beech et al. (2022) also show an increase in future simulated eddy activity that moves poleward in the ACC region; suggesting it to be driven by the intensification shifting westerly winds.

Fig. 9
figure 9

a Monthly SSH times series of the ACC region south of Australia for the historical runs. The two bold lines in black and purple show the timeseries of 2 random ensemble members from the 100-memew fibers. b Same as (a), but for the last 20 years from RCP8.5. c Monthly changes in standard deviation in meters for historical, RCP4.5 and RCP8.5 periods. A high pass filter 36 months has been applied to this timeseries

Full size image

To further investigate the variability changes and underlying mechanisms, we high pass filter the SSH variability with a 36-month window, focusing on the central location 106°E, 53°S. Results reveal that the intra-annual to interannual variability increases consistently from 12.5 cm for the historical to 15.0 and 19.5 cm for the RCP4.5 and RCP8.5 run, respectively (Fig. 9c). As it will be difficult to entail in detail the causes for the reduction and recovery of the low frequency variability here, we focus instead on the high frequency component as it can be linked directly to wind stress changes. Since the pattern of enhanced SSH variability SW of Australia shown in Fig. 4b has a similarity to the first mode of SSH variance shown by Webb and Cuevas (2002, their Fig. 1), their proposed relation to a trapped geostrophic mode serves as a starting point to also explain the change in variability seen here. According to those authors, wind driven Ekman pumping excites the first empirical mode, which they demonstrated through a correlation of up to 0.5 between area averaged wind stress curl and the first mode time series.

Fig. 10
figure 10

a Annual timeseries of easterly wind stress for the selected region in the Southern Ocean. The ensemble mean and the annual cycle are retained. b Monthly ensemble mean spread in easterly mean stress. A 60-month (5 year) running mean has been applied to this timeseries

Full size image

Following the close connection between wind stress curl and SSH, we perform a one-point correlation analysis across the ensemble members between the standard deviations of SSH and wind stress curl for RCP4.5 and RCP8.5 (Fig. 11). The data period is the last 20 years, and the time series were high pass (< 36 months) filtered to focus on the higher frequency part of the signal. The relation between curl and SSH is valid for the high frequency component. With the selected 36 month, enough samples are available to estimate variability sufficiently while the contribution from inter annual to decadal variability is small for which the relation increasingly holds less because of the increasing contribution from propagating signals. Note that the relevant signal on time scales shorter than 1 month is not resolved since only monthly mean data is available. For both scenarios, an enhanced correlation (> 0.2) is visible near and west of the location of the SSH variability marked by a green dot, confirming local Ekman pumping variability (shown as contours) as the relevant agent for modulating the SSH variability. This correlation is small due to analysis carried over a homogenous ensemble in which members have only small differences in statistics. When correlations are calculated from an ensemble containing both the historical and RCP8.5 members, then the correlations reach r = 0.85, which suggest that changes are predominantly driven by wind curl changes.

Also shown in Fig. 11 are the changes in the standard deviation of the Ekman pumping as difference between the RCPs minus the historical period. Both RCP runs show a southward shift of pumping variability SW of Australia signified by a zonally elongated dipole with the axis located roughly along 48°S. The shift leads to an enhanced pumping variability in a region that matches the enhanced SSH variability seen in Fig. 4b and provides therefore a simple mechanism for the enhanced variability. Furthermore, the change in Ekman pumping variability is about 2–3 times larger for RCP8.5 in comparison to RCP4.5 which agrees well with the change of high frequency SSH variability in this region.

We suggest the above SSH correlation in the southern hemisphere (Fig. 11) to be linked with changes in the Southern Ocean circulation and its influence on surface wind stress (and curl) that affects the Southern Ocean variability, through its southward movement of the westerly wind belt that circles Antarctica (see also Zheng et al. (2013). Most CMIP5 RCP scenarios suggest a trend in westerly winds towards stronger and more southward jets, correlated to a positive trend in the SAM (driven by opposing changes in stratospheric ozone concentrations and increasing concretions of GHGs) (Meijers 2014; Meijers et al. 2012).

Fig. 11
figure 11

One point correlation map showing correlations between standard deviations of SSH at 106°E, 53°S (marked by the green dot in both panels) and wind stress curl at every (global) grid point (high pass filter < 36 months) based on all month over 20 years periods. The black contours denote the change in wind stress curl standard deviation. Correlations are calculated across the ensemble dimension

Full size image

The MPI-GE wind stress changes for the high emission scenario of RCP8.5 agrees with the above findings in suggesting a projected pattern that resembles the positive phase of the SAM (Fig. 12). The threshold near the maximum RCP8.5 changes (see Figs. 9 and 10) for SSH variability and zonal wind stress is further substantiated by increased RCP8.5 mean zonal wind stress. It is 2.5 times larger compared to the RCP4.5 and shows characteristics of the positive phase of the SAM southwest of Australia (Fig. 12), suggesting that the wind variability shifts as the mean does. In summary, we suggest that this shift in the location of the wind is relevant, by the associated shift in wind stress curl variability. Here, we also suggest that future large scale atmospheric changes in the Southern Hemisphere and wind stress curl changes and respective ocean responses, induce this large RCP8.5 SSH threshold behavior. Several studies also showed increases in SSH variability linked to non-linear changes of threshold crossings seen from surface ocean warming (Widlansky et al. 2020). Identifying the exact mechanisms behind this is beyond the scope of this paper and requires an analysis of several single model large ensembles to detangle the mechanisms in the Southern Ocean responsible for such RCP8.5 SSH and zonal wind stress changes from mid-century onwards.

Fig. 12
figure 12

Difference in zonal wind stress between RCP 4.5-Hist and RCP 8.5-Hist for the SSH ensemble mean averaged over the 20 years periods 2081 to 2100 and 1986 to 2005, respectively. Stippling indicates significance at 95% level between the two periods based on a Student’s t-test

Full size image

5 Concluding remarks

In this paper we used the MPI-GE to investigate sterodynamic sea-level statistics and its changes under RCP4.5 and RCP8.5 climate change conditions. The large ensemble allows us to infer the sterodynamic sea-level statistics to an extent not previously conducted, thereby obtaining a more complete picture of sterodynamic sea-level statistics than hitherto available from climate models. Using the model output, we were able to investigate the nature of large-scale sea level variability represented by the model and its relation to Gaussian statistics.

The most pronounced deviation from Gaussian sterodynamic sea-level statistics exist in the model’s western and eastern tropical Pacific. However, by the end of 2100, SSH variability of the western tropical Pacific shows the tendency toward a more Gaussian variability under global warming conditions, linked in part with weaker zonal easterly wind stress pulses, suggesting a reduced El Niño Southern Oscillation activity in the western warm pool region in the future.

During the 21st century, regional changes in sterodynamic sea-level variability are most pronounced for the RCP 8.5 scenario, where sterodynamic sea-level variability is projected to increase in amplitude over much of the globe, with the largest increase in the ACC region visible southwest of Australia in the Southern Ocean correlated with changes in the wind stress curl which can lead to changes in the gyre circulation and shifts in the ACC. Similar distinct east-west changes are seen between SST variability and sterodynamic sea-level variability (Fig. 8) for the North Atlantic where an east (increase)-west (decrease) pattern emerges that can be influence from a combination of gyre strength (e.g., a weakened Gulf Stream), the declining AMOC strength and variability and the northeastward shifts in the westerly wind driven NAO and future evolution in ocean eddies. The combination of these complex mechanisms makes it challenging to attribute changes in sterodynamic sea-level variability to only one key mechanism. However, it is likely that the MPI-GE and its given resolution are associated with atmospherically forced variability and it changes under climate forcing.

Our results also suggest that extreme monthly sea level (computed as 95th and 99th percentiles) shift toward higher values than those expected from the pure change of the mean, which is linked primarily to changes in sea level variability. We diagnosed regional changes in the tails of the probability density functions (95th and 99th percentiles) with a global mean increase of 16 cm for RCP4.5 and by 24 cm for RCP8.5, respectively, suggesting increased high-end monthly sea level for warmer climate conditions.

Following from here, a caveat in our statistical sterodynamic sea-level analysis arise from uncertainties in the simulated sterodynamic sea-level changes. Since this model configuration is run at low resolution without resolving ocean eddies and smaller scale dynamics, not all processes affecting sea level are included in the model, such as storm surges and tides. This impacts the model’s response to forcing changes in dynamically active regions such as the Gulf Stream or the ACC southeast of Australia. Concerning eddies, it is vital to note that recent high-resolution simulations from Beech et al. (2022) predicts poleward intensification in eddy activity; albeit, with mixed results regionally. But one also needs to consider that the existing internal variability in the 100-member ensemble will be model specific and may be different when using another set of model ensembles. It therefore appears vital that our analysis gets repeated with different large ensemble results, but also with eddy-resolving models to either confirm that they are robust or identify their sensitivity to additional physics.

Even though large ensembles are limited to coarse resolution models, our study sheds light on implications on regional areas. The potential implications for coastal regions suggest increases in sterodynamic sea-level seen from the widening of the PDF tails, and increased variability along Western and North Pacific as well as Eastern North Atlantic. As we see sterodynamic sea-level variability increase in these regions, they will be prone to enhanced extreme sea level values due to internal ocean (eddy) variability that would be superimposed to a global mean SSH increase. If eddy variability would increase in the same way as we see the large variability to change here, natural variability can be a significant cause of enhanced high-end sterodynamic sea-level variability around the ocean. To verify this, in the future we need large ensembles based on eddy resolving climate models. It is vital to consider that many model ensemble outputs are needed to detect and attribute various statistics as being influenced by internal variability or distinguished as a forced response.