### Climate model simulations with CESM

In this study, we use model output of the Community Earth System Model (CESM) version 1.04. The high-resolution coupled climate version^{34} has an ocean component and sea-ice component with a (0.1^{circ }) horizontal resolution on a curvilinear, tri-polar grid which captures the development and interaction of mesoscale eddies^{11}. The ocean model (the Parallel Ocean Program) has 42 non-equidistant depth levels, with the highest vertical resolution near the surface. The atmosphere- and land surface components have a horizontal resolution of (0.5^{circ }), and the atmosphere component has 30 non-equidistant pressure levels.

We analyse two simulations of the high-resolution CESM which are initiated after a spin-up period of 200 years. The spin-up of the CESM has fixed forcing conditions ((hbox {CO}_2), methane, solar insolation, aerosols) of the year 2000 which are repeated every year. The two simulations consists of a control simulation (extension of the spin-up) and a forced simulation with a prescribed increase of atmospheric (hbox {CO}_2) levels (2000–2100) which is retained from the Representative Concentration Pathways (RCP) 8.5 scenario^{35}. The forced simulation is similar to a (1%) increase of (hbox {pCO}_2) each year. Both simulations have a 101 years’ duration (model years 200–300).

Supplementary Figure S1 shows some key variables for the spin-up period (similar results as in^{36}), the control simulation and the forced simulation for the high-resolution model. For the control simulation, the global mean (2 m) surface temperature is almost constant and the radiative imbalance is slightly positive. The Gregory plot (Supplementary Figure S1c) also shows the equilibration of the control simulation and the deviation from the equilibrium for the forced simulation. The upper 700 m global ocean heat content (Supplementary Figure S1d) is still adjusting, but relatively small trends (compared to the spin-up) occur over the last 100 years. The deep ocean fields take a much longer time to equilibrate. In the forced simulation, the surface temperature, radiative imbalance and ocean heat content start to deviate from the control simulation due to the increase of atmospheric (hbox {CO}_2) concentration.

In addition to the high-resolution CESM simulations, a companion simulation was conducted at a lower resolution (using CESM version 1.1.2). This model’s ocean component has a horizontal resolution of (1^{circ }) and 60 non-equidistant depth levels. The atmospheric component has a horizontal resolution of (1^{circ }) and 30 non-equidistant pressure levels. The low-resolution simulation of the CESM had a spin-up period of 1,200 years with the same fixed forcing conditions as the high-resolution version. After the spin-up, we conduct the two simulations: one control simulation and one forced simulation, similar to the set-up of the high-resolution version. Both simulations have a 101 years’ duration (model years 1,200–1,300).

### Model output analysed

The standard model output of the CESM simulations consists of monthly averaged fields of sea surface height above geoid ((eta _M)), horizontal velocity, temperature and salinity. All of the monthly averaged quantities are converted to yearly averages. The globally averaged (eta _M) fields are about zero since the ocean is volume conserving. A limited number of quantities, for example (eta _M), are available as daily averages for the CESM simulations and are used in generalised extreme value analysis.

### CMIP6 data

We use results from the latest release of the Coupled Model Intercomparison Project phase six (CMIP6) and compare these to the output of our CESM simulations (Table 1). We analysed the preliminary model output of the CMIP6 in which atmospheric (hbox {CO}_2) levels increase each year by (1%). Note that each model of the CMIP6 is initiated from the pre-industrial (year 1,850) control simulation. Only models which conserve ocean volume (as the CESM) are considered here (variable name ‘zos’). We analyse the monthly averaged (eta _M), horizontal velocity, temperature and salinity fields of the first 101 model years, as is done for the CESM output. The model output of the AWI-CM-1-1-MR is provided on an unstructured grid. Before analysing the output of this CMIP6 model, the data is interpolated onto the HR-CESM grid.

### Significance of the linear trends

The trends are derived from a linear least-square fit to the yearly-averaged time series. The significance of each trend is determined following the procedure outlined in^{37}, while taking into account the reduction of degrees of freedom for the time series which are not statistically independent. First, the variance (s_e) of the residuals of the linear fit (*e*(*t*)) are determined. The degrees of freedom are reduced using the lag-1 autocorrelation ((r_1)) of the residuals,

$$begin{aligned} s_e^2 = frac{1}{N frac{1 – r_1}{1 + r_1} – 2} sum _{t = 1}^N e(t)^2 end{aligned}$$

(1)

where *N* is the number of years (i.e. (N = 101) years). The standard error of the residuals, (s_b) is

$$begin{aligned} s_b = frac{s_e}{sqrt{sum _{t = 1}^{N} left( t – frac{N + 1}{2} right) ^2 }} end{aligned}$$

(2)

The Student-*t* value is the ratio between the linearly fitted trend and the standard error. Using the reduced degrees of freedom and the two-sided critical Student-*t* values, one can determine the significance of having a trend different from zero (the null hypothesis).

### Sterodynamic sea level

The local sterodynamic sea level (SDSL) consists of two components^{22}. The first component is the dynamic sea level ((eta _M)) and is the height of the ocean surface above the geoid. The (eta _M) fields are part of the standard output in the CESM simulations (variable ‘SSH’) and CMIP6 models (variable ‘zos’). The (eta _M) fields have a global mean of about 0, if the global mean was non zero we subtracted uniformly the global mean from the (eta _M) fields.

The next component is the global-mean thermosteric sea-level rise ((eta _S^g)), which is determined from post-processing the model output^{38}. First, we determined the local steric sea-level ((eta _S)). The contribution of both thermal changes and haline changes is determined as the full-depth integral over the specific volume anomaly^{39},

$$begin{aligned} eta _S = int _{-H}^{0} frac{rho _0 – rho (T, S, P)}{rho _0} mathrm {d}z end{aligned}$$

(3)

The temperature, salinity and pressure dependency are taken into account while determining the density, and (rho _0 = 1,028) kg m(^{-3}). The steric sea-level change is expressed as an anomaly with respect to the initial value of the first model year. The global-mean thermosteric sea-level ((eta _S^g)) rise is determined by globally averaging (eta _S). This procedure is done for both the CESM simulations as the CMIP6 models. One can use the variable ‘zostoga’ instead of (eta _S^g) for the CMIP6 models, but this is not used in this study for comparison with the CESM simulations.

Eventually, hence the local sterodynamic sea level ((eta)) is determined from

$$begin{aligned} eta = eta _M + eta _S^g end{aligned}$$

(4)

The trend over the 101-year period for (eta _M) and (eta _S) are shown in Supplementary Figure S3. We corrected for any drift in (eta _S) and (eta _S^g) using the control simulations.

### Barotropic streamfunction

The barotropic flow is defined as the full depth integral of the horizontal velocity:

$$begin{aligned} {overrightarrow{BF}} = int _{-H}^{0} vec {v},mathrm {d}z end{aligned}$$

(5)

Starting from Antarctica (with a value of 0 for the barotropic streamfunction), we integrate the zonal component of the barotropic flow (indicated by (BF_x)) meridionally to determine the barotropic streamfunction:

$$begin{aligned} mathrm {BSF} = int _{90^{circ }mathrm {S}}^{90^{circ }mathrm {N}} BF_x,mathrm {d}y end{aligned}$$

(6)

For convenience, the average value of the BSF along the African coast line is subtracted from the entire BSF field.

### AMOC strength (definition)

The AMOC strength is defined as the total meridional mass transport at (26^{circ },hbox {N}) integrated between (80.5^{circ },hbox {W}) and (12^{circ },hbox {W}) (RAPID array) and integrated over the upper 1,000 m:

$$begin{aligned} mathrm {AMOC} = int _{-1000}^{0} int _{80.5^{circ }mathrm {W}}^{12^{circ }mathrm {W}} v,mathrm {d}x mathrm {d}z end{aligned}$$

(7)

Observations show that the AMOC strength (3-month average) at (26^{circ },hbox {N}) varies between 10 and 24 Sv (1 Sv (equiv)(10^6,hbox {m}^3,hbox {s}^{-1})), with a mean transport of about 16–19 Sv^{21,40}.

### NADW strength (definition)

The NADW strength is defined as the total meridional mass transport at (40^{circ },hbox {N}) integrated between (77^{circ },hbox {W}) and (8^{circ },hbox {W}) and between 1500 and 4000 m.

$$begin{aligned} mathrm {NADW} = int _{-4000}^{-1500} int _{77^{circ }mathrm {W}}^{8^{circ }mathrm {W}} v,mathrm {d}x mathrm {d}z end{aligned}$$

(8)

### Generalised extreme value

The (eta _M) extremes in the NBC outflow region are analysed using GEV (Generalised Extreme Value) analysis. First, we retained the monthly (eta _M) maxima ((eta _M^{Max})) from daily averaged (eta _M). This is the block maxima GEV approach where each month contains the (eta _M) maximum. The distribution of (eta _M^{Max}) may vary in time. Therefore different periods of *n* months’ length are retained from the full time series, under the assumption that these periods in time are stationary and stationary GEV analysis can be conducted^{41}. Using all the data points over a particular period, a GEV fit is made to the data points using the following expression as the cumulative distribution function:

$$begin{aligned} G_n(z) = exp left( – left( 1 + xi left( frac{z – mu }{sigma } right) right) ^{-1 / xi } right) end{aligned}$$

(9)

where (mu), (sigma), (xi) are the location, scale and shape parameter, respectively. A special case of the distribution is where (xi rightarrow 0), the GEV distribution approaches a Gumbel distribution. The probability of occurrence (*p*) using the distribution of (G_n(z)) is related to return time ((tau)) of such an event:

$$begin{aligned} tau = frac{1}{p} end{aligned}$$

(10)

The value of (eta _M^{Max}) related to this probability ((z_p)) can be obtained by:

$$begin{aligned} z_p = mu – frac{sigma }{xi } left( 1 – left( – log left( 1 – p right) right) ^{-xi } right) end{aligned}$$

(11)

Note that detrending the periods of time (linear or quadratic) might result in biases. Therefore, the data on which GEV is applied is directly fitted to the (G_n(z)) distribution. We have chosen to fit the (G_n(z)) distribution using periods of 360 months’ length (30 years). Shorter sections (120 months) provided similar results, but might have insufficient data to make a reasonable fit.