Warning of a forthcoming collapse of the Atlantic meridional overturning circulation
The Atlantic meridional overturning circulation (AMOC) is a major tipping element in the climate system and a future collapse would have severe impacts on the climate in the North Atlantic region. In recent years weakening in circulation has been reported, but assessments by the Intergovernmental Panel on Climate Change (IPCC), based on the Climate Model Intercomparison Project (CMIP) model simulations suggest that a full collapse is unlikely within the 21st century. Tipping to an undesired state in the climate is, however, a growing concern with increasing greenhouse gas concentrations. Predictions based on observations rely on detecting early-warning signals, primarily an increase in variance (loss of resilience) and increased autocorrelation (critical slowing down), which have recently been reported for the AMOC. Here we provide statistical significance and data-driven estimators for the time of tipping. We estimate a collapse of the AMOC to occur around mid-century under the current scenario of future emissions. The Atlantic meridional overturning circulation (AMOC) is a major tipping element in the climate system. Here, data-driven estimators for the time of tipping predict a potential AMOC collapse mid-century under the current emission scenario.
/static/images/favicons/nature/apple-touch-icon-f39cb19454.png
Nature
A forthcoming collapse of the Atlantic meridional overturning circulation (AMOC) is a major concern as it is one of the most important tipping elements in Earth’s climate system1,2,3. In recent years, model studies and paleoclimatic reconstructions indicate that the strongest abrupt climate fluctuations, the Dansgaard-Oeschger events4, are connected to the bimodal nature of the AMOC5,6. Numerous climate model studies show a hysteresis behavior, where changing a control parameter, typically the freshwater input into the Northern Atlantic, makes the AMOC bifurcate through a set of co-dimension one saddle-node bifurcations7,8,9. State-of-the-art Earth-system models can reproduce such a scenario, but the inter-model spread is large and the critical threshold is poorly constrained10,11. Based on the CMIP5 generation of models, the AR6 IPCC report quotes a collapse in the 21st century to be very unlikely (medium confidence)12. Among CMIP6 models, there is a larger spread in the AMOC response to warming scenarios, thus an increased uncertainty in the assessment of a future collapse13. There are, however, model biases toward overestimated stability of the AMOC, both from tuning to the historic climate record14, poor representation of the deep water formation15, salinity and glacial runoff16.

When complex systems, such as the overturning circulation, undergo critical transitions by changing a control parameter λ through a critical value λc, a structural change in the dynamics happens. The previously statistically stable state ceases to exist and the system moves to a different statistically stable state. The system undergoes a bifurcation, which for λ sufficiently close to λc can happen in a limited number of ways rather independent from the details in the governing dynamics17. Besides a decline of the AMOC before the critical transition, there are early-warning signals (EWSs), statistical quantities, which also change before the tipping happens. These are critical slowing down (increased autocorrelation) and, from the Fluctuation-Dissipation Theorem, increased variance in the signal18,19,20. The latter is also termed “loss of resilience”, especially in the context of ecological collapse21. The two EWSs are statistical equilibrium concepts. Thus, using them as actual predictors of a forthcoming transition relies on the assumption of quasi-stationary dynamics.

The AMOC has only been monitored continuously since 2004 through combined measurements from moored instruments, induced electrical currents in submarine cables and satellite surface measurements22. Over the period 2004–2012, a decline in the AMOC has been observed, but longer records are necessary to assess the significance. For that, careful fingerprinting techniques have been applied to longer records of sea surface temperature (SST), which, backed by a survey of a large ensemble of climate model simulations, have found the SST in the Subpolar gyre (SG) region of the North Atlantic (area marked with a black contour in Fig. 1a) to contain an optimal fingerprint of the strength of the AMOC23,24,25.

Figure 1b shows the SG and the GM SSTs obtained from the Hadley Centre Sea Ice and Sea Surface Temperature data set (HadISST)26. Figure 1c shows the SG anomaly, and Fig. 1d shows the GM anomaly with a clear global warming trend in the last half of the record. The AMOC fingerprint for the period 1870–2020 is shown in Fig. 1e. This is the basis for the analysis. It has been reported11,27 that this and similar AMOC indices show significant trends in the mean, the variance and the autocorrelation, indicating early warning of a shutdown of the AMOC. However, a trend in the EWSs within a limited period of observation could be a random fluctuation within steady-state statistics. Thus, for a robust assessment of the shutdown, it is necessary to establish a statistical confidence level for the change above the natural fluctuations. This is not easily done given only one, the observed realization of the approach to the transition. Here we establish such a measure of the confidence for the variance and autocorrelation and demonstrate that variance is the more reliable of the two. A further contribution is an estimator of not only whether a transition is approaching but also the time when the critical transition is expected to occur. The strategy is to infer the evolution of the AMOC solely on observed changes in mean, variance and autocorrelation. The typical choice of control parameter is the flux of freshwater into the North Atlantic. River runoff, Greenland ice melt and export from the Arctic Ocean are not well constrained28; thus, we do not assume the control parameter known. Boers27 assumes the global mean temperature T to represent the control parameter. Although T has increased since ~1920 (Fig. 1d), the increase is not quite linear with time. All we assume here is that the AMOC is in an equilibrium state prior to a change toward the transition. The simplest uninformed assumption is that the change is sufficiently slow and that the control parameter approaches the (unknown) critical value linearly with time. This assumption is confirmed by a close fit of the estimated model to the observed AMOC fingerprint. Although we make no explicit assumptions, the primary driver of climate change, the logarithm of the atmospheric CO2 concentration, does, in fact, increase close to linearly with time in the industrial period29. Our results are robust without making specific assumptions regarding the driver of the AMOC.

In this work, we show that a transition of the AMOC is most likely to occur around 2025-2095 (95% confidence interval).

Modeling and detecting the critical transition

Denote the observed AMOC fingerprint by x(t) (Fig. 1e). We model it by a stochastic process Xt, which, depending on a control parameter λ < 0, is at risk of undergoing a critical transition through a saddle-node bifurcation for λ = λc = 0. The system is initially in a statistically stable state, i.e., it follows some stationary distribution with constant λ = λ0. We are uninformed about the dynamics governing the evolution of Xt but can assume effective dynamics, which, with λ sufficiently close to the critical value λc = 0, can be described by the stochastic differential equation (SDE):

$$d{X}_{t}=-(A{({X}_{t}-m)}^{2}+\lambda )dt+\sigma d{B}_{t},$$

(1)

where \(m=\mu -\sqrt{|\lambda|/A}\) and μ is the stable fixed point of the drift, A is a time scale parameter, Bt is a Brownian motion and σ2 scales the variance. Disregarding the noise, this is the normal form of the co-dimension one saddle-node bifurcation17 (see “Methods”). The square-root dependence of the stable state: \(\mu -m \sim \sqrt{{\lambda }_{c}-\lambda }\) is the main signature of a saddle-node bifurcation. It is observed for the AMOC shutdown in ocean-only models as well as in coupled models, see Fig. 2, in strong support of Eq. (1) for the AMOC.

At time t0, λ(t) begins to change linearly toward λc = λ(tc) = 0:

$$\lambda (t)={\lambda }_{0}(1-{{\Theta }}[t-{t}_{0}](t-{t}_{0})/{\tau }_{r}),$$

(2)

where Θ[t] is the Heaviside function and τr = tc − t0 > 0 is the ramping time up to time tc, where the transition eventually will occur. Time tc is denoted the tipping time; however, an actual tipping can happen earlier than tc due to a noise-induced tipping. As the transition is approached, the risk of noise-induced tipping (n-tipping) prior to tc is increasing and, at some point, making the EWSs irrelevant for predicting the tipping. The probability for n-tipping can, in the small noise limit, be calculated in closed form, \(P(t,\lambda )=1-\exp (-t/{\tau }_{n}(\lambda ))\), with mean waiting time \({\tau }_{n}(\lambda )=(\pi /\sqrt{|\lambda|})\exp (8|\lambda {|}^{3/2}/3{\sigma }^{2})\) (see “Methods”).

The mean and variance are calculated from the observations as the control parameter λ(t) is possibly changing. These EWSs are inherently equilibrium concepts and statistical; thus, a time window, Tw, of a certain size is required for a reliable estimate. As the transition is approached, the differences between the EWSs and the pre-ramping values of the variance and autocorrelation (baseline) increase; thus, a shorter window Tw is required for detecting a difference. Conversely, close to the transition critical slowing down decreases the number of independent points within a window, thus calling for a larger window for reliable detection. Within a short enough window, [t − Tw/2, t + Tw/2], we may assume λ(t) to be constant and the noise small enough so that the process (1) for given λ is well approximated by a linear SDE, the Ornstein–Uhlenbeck process30. A Taylor expansion around the fixed point μ(λ) yields the approximation

$$d{X}_{t}\approx -\alpha (\lambda )({X}_{t}-\mu (\lambda ))dt+\sigma d{B}_{t}$$

(3)

where \(\mu (\lambda )=m+\sqrt{|\lambda|/A}\) and \(\alpha (\lambda )=2\sqrt{A|\lambda|}\) is the inverse correlation time. For fixed λ, the process is stationary, with mean μ, variance γ2 = σ2/2α and one-lag autocorrelation \(\rho=\exp (-\alpha {{\Delta }}t)\) with step size Δt = 1 month. As λ(t) increases, α decreases, and thus variance and autocorrelation increase. From μ, γ2 and ρ the parameters of Eq. (1) are determined: \(\alpha=-\log \rho /{{\Delta }}t\), σ2 = 2αγ2, A = α/2(μ − m) and \(\lambda={({\sigma }^{2}/4{\gamma }^{2})}^{2}/A\). Closed form estimators for μ, γ2 and ρ are obtained from the observed time series within a running window by maximum likelihood estimation (MLE) (Supplementary text S1, see also ref. 31).

The uncertainty is expressed through the variances of the estimators \({\hat{\gamma }}^{2}\) and \(\hat{\rho }\) obtained from the observations within a time window Tw. The hats indicate that they are estimators and thus stochastic variables with variances around the true values. Detection of an EWS at some chosen confidence level q (such as 95 or 99%) requires one of the estimates \({\hat{\gamma }}^{2}\) or \(\hat{\rho }\) for a given window to be statistically different from the baseline values \({\hat{\gamma }}_{0}^{2}\) or \({\hat{\rho }}_{0}\), which depend on the window size as well as how different the EWSs are from their baseline values.

Time scales in early-warning signals

The detection of a forthcoming transition using statistical measures involves several time scales. The primary internal time scale is the autocorrelation time, tac, in the steady state. The ramping time τr over which the control parameter changes from the steady state value to the critical value sets an external time scale. For given α(λ) and q-percentile, the required time window Tw(q, α) to detect a change from baseline in EWSs at the given confidence level q is given in the closed form in the next section (Eq. (7) for variance and Eq. (8) for autocorrelation). The approach to the collapse and the involved time scales are schematically summarized in Fig. 3, while they are calculated in Fig. 4a, where the required window size Tw at the 95% confidence level is plotted as a function of λ for the variance (red curve) and autocorrelation (yellow curve). These are plotted together with the mean waiting time for n-tipping, tnoise, (blue curve). With Tw = 50 years, increased variance can only be detected after the time when λ(t) ≈ −1.2 (crossing of red and red-dashed curves). At that time, a window of approximately 75 years is required to detect an increase in autocorrelation, making variance the better EWS of the two. When λ ≈ −0.4, the mean waiting time for n-tipping is smaller than the data window size. Thus, the increased variance can be used as a reliable EWS in the range −1.2 < λ(t) < −0.4, indicated by the green band. How timely an early warning this is depends on the speed at which λ(t) is changing from λ0 to λc, i.e., the ramping time τr. A set of 1000 realizations has been simulated with λ0 = −2.82 and τr = 140 years, indicated by the time labels on top of Fig. 4a. Ten of these realizations are shown in Fig. 4b on top of the stable and unstable branches of fixed points of model (1) (the bifurcation diagram). Figure 4c (d) shows the variance (autocorrelation) calculated from the realizations within a running 50-year window (shown in Fig. 4c). The solid black line is the baseline value for λ = λ0, while the solid blue line is the increasing value for λ = λ(t). The calculated 95% confidence level for the measurement of the EWS within the running window is shown by the dashed black and blue lines, respectively. The corresponding light blue curves are obtained numerically from the 1000 realizations. The green band in Fig. 4c corresponds to the green band in Fig. 4a and shows where early warning is possible in this case.

Statistics of early-warning signals

The variances of the estimators are approximately (see Supplementary text S1).

$${{{{{{{\rm{Var}}}}}}}}({\hat{\gamma }}^{2})\, \approx \, \frac{2{({\gamma }^{2})}^{2}}{\alpha {T}_{w}}=\frac{{\sigma }^{4}}{2{\alpha }^{3}{T}_{w}};\,{{{{{{{\rm{Var}}}}}}}}(\hat{\rho })\, \approx \, \frac{2\alpha {{\Delta }}{t}^{2}}{{T}_{w}},$$

(4)

where Tw = nΔt is the observation window.

The question is then how large Tw needs to be to detect a statistically significant increase compared to the estimated baseline values \({\hat{\gamma }}_{0}^{2}\) and \({\hat{\rho }}_{0}\). For a given estimate \({\hat{\gamma }}^{2}\), the estimated difference from the baseline variance is

$${{{\Delta }}}_{{\gamma }^{2}}={\hat{\gamma }}^{2}-{\hat{\gamma }}_{0}^{2}={\hat{\gamma }}_{0}^{2}({\hat{\alpha }}_{0}/\hat{\alpha }-1),$$

(5)

and the estimated difference from the baseline autocorrelation is

$${{{\Delta }}}_{\rho }=\hat{\rho }-{\hat{\rho }}_{0}={\hat{\rho }}_{0}({e}^{({\hat{\alpha }}_{0}-\hat{\alpha }){{\Delta }}t}-1)\, \approx \, {\hat{\rho }}_{0}({\hat{\alpha }}_{0}-\hat{\alpha }){{\Delta }}t.$$

(6)

Since the two EWSs, \({\hat{\gamma }}^{2}\) and \(\hat{\rho }\), are treated on an equal footing, in the following, we let \(\hat{\psi }\) denote either of the estimators (given explicitly in Supplementary text S1, Eqs. (S5) or (S6)). The standard error is \(s(\hat{\psi })={{{{{{{\rm{Var}}}}}}}}{(\hat{\psi })}^{1/2}\) (Eq. (4)) and \(\hat{{{\Delta }}}\) denotes either of the two estimated differences (5) or (6). The null hypothesis is that λ = λ0, or equivalently α = α0. The null distribution of \(\hat{\psi }\) is assumed to be Gaussian (confirmed by simulations). A quantile q from the standard Gaussian distribution expresses the acceptable uncertainty in measuring the statistical quantity ψ. We thus get that \(\hat{{{\Delta }}} \, < \, qs(\hat{\psi })\) at the q-confidence level (95%, 99% or similar) under the null hypothesis. To detect an EWS at the q-confidence level based on measuring ψ at time t, we require that \(\hat{{{\Delta }}}(t) \, > \, q(s(\hat{\psi }(t))+s({\hat{\psi }}_{0}))\), which, solved for Tw gives for variance:

$${T}_{w} > 2{q}^{2}{\left(\frac{\hat{\alpha }(t)/\sqrt{{\hat{\alpha }}_{0}}+{\hat{\alpha }}_{0}/\sqrt{\hat{\alpha }(t)}}{{\hat{\alpha }}_{0}-\hat{\alpha }(t)}\right)}^{2},$$

(7)

and for autocorrelation,

$${T}_{w} > 2{q}^{2}{\left(\frac{\sqrt{{\hat{\alpha }}_{0}}+\sqrt{\hat{\alpha }(t)}}{{\hat{\alpha }}_{0}-\hat{\alpha }(t)}\right)}^{2}{\hat{\rho }}_{0}^{-2}.$$

(8)

Substituting \({\alpha }_{0}=2\sqrt{A|{\lambda }_{0}|}\) and \(\alpha (t)=2\sqrt{A|\lambda (t)|}\) provides the time window Tw needed to detect an EWS at time t with large probability. Eqs. (7) and (8) are illustrated in Fig. 4a (red and yellow curves), where it is seen that detecting a significant increase in variance requires a shorter data window than detecting a significant increase in autocorrelation. Two times \(s(\hat{\psi }(t))\) around the mean of the ramped variance and two times \(s({\hat{\psi }}_{0})\) around baseline values are illustrated in Fig. 4c, d (dashed lines). Once a trace leaves the baseline confidence interval, a statistically significant change is detected, and when the two dashed lines cross, 95% of the traces have detected an EWS (Eqs. (7) and (8)).

Predicting a forthcoming collapse of the AMOC

The AMOC fingerprint shown in Fig. 1e (replotted in Fig. 5a) shows an increased variance, γ2, and autocorrelation, ρ, plotted in Fig. 5b, c as functions of the mid-point of a 50-year running window, i.e., the EWS obtained in 2020 is assigned to the year 1995. The estimates leave the confidence band of the baseline values (pink area) around the year 1970. This is not the estimate of t0, which happened earlier and is still to be estimated; it is the year where EWSs are statistically different from baseline values. The estimates after 1970 stay consistently above the upper limit of the confidence interval and show an increasing trend, and we thus conclude that the system is moving toward the tipping point with high probability.

To estimate the tipping time once it has been established that the variance and autocorrelation are increasing, we use two independent methods to check the robustness of our results: (1) Moment-based estimator that uses the variance and autocorrelation estimates within the running windows. (2) Approximate MLE directly on model (1)-(2) with no running window. The advantage of the first method is that it has less model assumptions; however, it is sensitive to the choice of window size. The advantage of the second method is that it uses the information in the data more efficiently given model (1)-(2) is approximately correct, it has no need for a running window and does not assume stationarity after time t0. In general, MLE is statistically the preferred method of choice, giving the most accurate results with the lowest estimation variance.

The first method, the moment estimator of the tipping time obtains, within the running window, the parameters α(t) (Fig. 5d) and σ2 (Fig. 5e) of the linearized dynamics, Eq. (3), and thus also γ2(t). Within the running window, the data are detrended before estimation by subtracting a linear regression fit in order not to falsely inflate the variance estimates caused by deviations from stationarity. Then we obtain Aλ(t) from σ2 and γ2(t) (Fig. 5f) using that \(A\lambda (t)={({\sigma }^{2}/4{\gamma }^{2}(t))}^{2}\). This is consistent with a linear ramping of λ(t) beginning from a constant level λ0 at a time t0. By sweeping t0 from 1910 to 1950 and Tw from 45 to 65 years, we obtain Aλ0 and τr from the least square error fit to the data. This shows a single minimum at t0 = 1924 and Tw = 55 years (Fig. 6e). Setting t0 = 1924, we obtain tc from a linear fit (regressing λ on t) from the crossing of the x-axis (λc = 0). This is shown in Fig. 5f (red line). This yields −Aλ0 = 2.34 year−2 and τr = 133 years. Thus, the tipping time is estimated to be in the year 2057, shown in Fig. 5f. Since we have only obtained the combined quantity \(A\lambda={({\sigma }^{2}/4{\gamma }^{2})}^{2}\), we still need to determine A and m in Eq. (1). We do that from the best linear fit to the mean level \(\mu=m+\sqrt{|\lambda|/A}\) observing that \(\mu=m+\sqrt{A|\lambda|}(1/A)=m+({\sigma }^{2}/4{\gamma }^{2})(1/A)\). The estimates are shown by the red curves in Fig. 5a–f. The red dot in Fig. 5a is the tipping point, and the dashed line in Fig. 5b is the asymptote for the variance. With the parameter values completely determined, the confidence levels are calculated: The two-standard error levels around the baseline values of the EWS are shown by purple bands in Fig. 5b, c. Thus, both EWSs show increases beyond the two-standard error level from 1970 and onward.

The second method, the approximate MLE of the tipping time, is applied to model (1)–(2). The likelihood function is the product of transition densities between consecutive observations. However, the likelihood is not explicitly known for this model, and we therefore approximate the transition densities. From the data before time t0, approximation (3) is used, where exact MLEs are available (Supplementary text S1). This provides estimates of the parameters λ0, m as a function of parameter A, as well as the variance parameter σ2. To estimate A and τr, the observations after time t0 are used. After time t0, the linear approximation (3) is no longer valid because the dynamics are approaching the bifurcation point, and the non-linear dynamics will be increasingly dominating. The likelihood function is the product of transition densities, which we approximate with a numerical scheme, the Strang splitting, which has shown to have desirable statistical properties for highly non-linear models, where other schemes, such as the Euler–Maruyama approximation is too inaccurate32 (Supplementary text S2). Using t0 = 1924, the optimal fit is the same as the moment method, tc = 2057, with a 95% confidence interval 2025–2095.

Confidence intervals for the estimate of the tipping time are obtained by bootstrap. The likelihood approach provides asymptotic confidence intervals; however, these assume that the likelihood is the true likelihood. To incorporate also the uncertainty due to the data generating mechanism (1) not being equal to the Ornstein–Uhlenbeck process (3) used in the likelihood, we chose to construct parametric bootstrap confidence intervals. This was obtained by simulating 1000 trajectories from the original model with the estimated parameters and repeating the estimation procedure on each data set. Empirical confidence intervals were then extracted from the 1000 parameter estimates. These were indeed larger than the asymptotic confidence intervals provided by the likelihood approach, however, not by much. Histograms of the bootstrapped estimates are shown in Fig. 6a–d. The histogram in Fig. 6a is the tipping year, repeated in yellow in 5f.

The mean of the bootstrapped estimates of the tipping time is 〈tc〉 = 2050, and the 95% confidence interval is 2025–2095. The small discrepancy in the mean is probably due to the approximate model used for estimation being different from the data-generating model (1), confirming that the linear model still provides valid estimates even if the true dynamics are unknown. To test the goodness-of-fit, normal residuals (see “Methods”) were calculated for the data. These are plotted in Fig. 6f as a quantile-quantile plot. If the model is correct, the points fall close to a straight line. The model is seen to fit the data well, further supporting the obtained estimates.

To obtain the AMOC fingerprint, two steps are required: The seasonal cycle in the SST is governed by the surface radiation independent from the circulation and thus removed by considering the monthly anomalies, where the mean over the period of recording of the month is removed. Second, there is an ongoing positive linear trend in the SST related to global warming, which is also not related to circulation. This is compensated for by subtracting 2 × the global mean (GM) SST anomaly (small seasonal cycle removed). This differs slightly from ref. 23, where 1 × the GM SST was subtracted. The “translation” from the proxy SST temperature and AMOC flow is 0.26 SV/K [ref. 23, Fig. 3]. Here we have taken into account that the warming is not globally homogeneous: The warming in the SG region is larger than the global mean due to polar amplification. The way we have estimated this effect is by comparing the proxy with the AMOC estimates covering the period 1957–2004 from the so-called MOCz reported in the review by ref. 42. This shows a drop of 3 SV in that period. Minimizing the difference between the proxy SSTSG-β SSTGM and this more direct measurement with respect to β, we get β = 1.95 ≈ 2 rather than β = 1 used by ref. 23. The factor 2 is thus the optimal value for the polar amplification43 obtained by calibrating to recent direct measurements42. The original and our calibrated proxies are shown in Fig. 7.

To check the robustness with respect to the AMOC fingerprint record, we repeated the analysis, subtracting 1x and 3x GM SST from the SG SST. Subtracting 3x GM SST only changes estimate and the confidence intervals by a few years, whereas subtracting 1x GM SST delays the tipping with 25 years, but the overall trend and conclusions do not change. The results are given in Table 1. In the reanalysis, we fixed t0 = 1924.

Estimator of the tipping time and model control

The AMOC fingerprint is assumed to be observations from a process Xt given as a solution to Eqs. (1) and (2), and we wish to estimate the parameters θ = (A, m, λ0, τr, σ) from observations (x0, x1, …, xn) before time t0 and observations (y0, y1, …, yn) after time t0. These equations cannot be explicitly solved, and the exact distribution of Xt is not explicitly known. A standard way to solve this is to approximate the transition density by a Gaussian distribution obtained by the Euler–Maruyama scheme. However, the estimators obtained from the Euler–Maruyama pseudo-likelihood are known to be biased, especially in non-linear models32. We estimate by a two-step procedure, approximating the stationary distribution before time t0 by an Ornstein–Uhlenbeck process, of which exact maximum likelihood estimators are available (see Supplementary text S1), and using Strang splitting for the non-stationary and non-linear part after time t0, using methods proposed in ref. 32, see Supplementary text S2 for details.

To test the model fit, uniform residuals, ui, i = 1, …, n were calculated for the AMOC data using the estimated parameters from the MLE method as follows. The model assumes that observation xi follows some distribution function \({F}_{i,\hat{\theta }}\) for the estimated parameter values \(\hat{\theta }\). If this is true, then \({u}_{i}={F}_{i,\hat{\theta }}({x}_{i})\) is uniformly distributed on (0, 1). Transforming these residuals back to a standard normal distribution provides standard normally distributed residuals if the model is true. Thus, a normal quantile-quantile plot reveals the model fit. The points should fall close to a straight line. The reason for making the detour around the uniform residuals is twofold. First, since the data is not stationary, each observation follows its own distribution, and residuals cannot be directly combined. Second, since the model is stochastic, standard residuals are not well-defined, and observations should be evaluated according to their entire distribution, not only the distance to the mean.

Noise-induced tipping

The drift term in Eq. (1) is the negative gradient of a potential, f(x, λ) = − ∂xV(x, λ) = − (A(x−m)2 + λ) with V(x, λ) = A(x−m)3/3 + (x − m)λ. For λ < 0, the drift has two fixed points, \(m\pm \sqrt{|\lambda|/A}\). The point \(m+\sqrt{|\lambda|/A}\) is a local minimum of the potential V(x, λ) and is stable, whereas \(m-\sqrt{|\lambda|/A}\) is a local maximum and unstable. The system thus has two basins of attraction separated by \(m-\sqrt{|\lambda|/A}\), with a drift toward either \(m+\sqrt{|\lambda|/A}\) or −∞ dependent on whether \({X}_{t} > m-\sqrt{|\lambda|/A}\) or \({X}_{t} < m-\sqrt{|\lambda|/A}\). We denote the two basins of attraction, the normal and the tipped state, respectively. When λ = 0, the normal state disappears, and the system undergoes a bifurcation and Xt will be drawn toward −∞.

Due to the noise, the process can escape into the tipped state by crossing over the potential barrier \({{\Delta }}(\lambda )=V(m-\sqrt{|\lambda|/A},\,\lambda )\,-\,V(m+\sqrt{|\lambda|/A},\,\lambda )=4|\lambda {|}^{3/2}/3{A}^{1/2}\). Assume Xt to be close to \(m+\sqrt{|\lambda|/A}\) at some time t, i.e., in the normal state. The escape time will asymptotically (for σ → 0) follow an exponential distribution such that

$$P(t,\lambda )=1-\exp (-t/{\tau }_{n}(\lambda ))$$

(9)

where P(t, λ) is the probability of observing an escape time shorter than t for a given value of λ. The mean noise-induced escape time τn(λ) is44,45:

$${\tau }_{n}(\lambda )=\frac{2\pi \exp (2{{\Delta }}(\lambda )/{\sigma }^{2})}{\sqrt{{V}^{{\prime\prime} }(m+\sqrt{|\lambda|/A},\lambda )|{V}^{{\prime\prime} }(m-\sqrt{|\lambda|/A},\lambda )|}} \\=(\pi /\sqrt{A|\lambda|})\exp (8|\lambda {|}^{3/2}/3{A}^{1/2}{\sigma }^{2}).$$

(10)

Assume that the rate of change of λ(t) follows Eq. (2), then for τr < τn(λ), the waiting time for a random crossing is so long that a crossing will not happen before a bifurcation-induced transition happens (b-tipping). If τr > τn(λ), a noise-induced tipping is expected before the bifurcation point is reached. Since τn(λ) decreases with increasing λ, at some point, the two time scales will end up matching.

Normal form of the saddle-node bifurcation

Consider the general dynamical equation

$$\frac{dx}{dt}=f(x,\lambda ),$$

(11)

where x is a variable and λ is a (fixed) parameter. A point x0 with f(x0, λ) = 0 is a fixed point or steady state. A fixed point is stable/unstable if \({\partial }_{x}f{(x,\lambda )}_{x={x}_{0}}\) is negative/positive; thus, the fixed point is attracting/repelling. If f(x, λ) is not a linear function of x, multiple steady states may exist. A saddle-node bifurcation occurs when changing the control parameter λ through a critical value λc a stable and an unstable fixed point merge and disappear. The situation is shown in the figure, where the blue surface is f(x, λ), while the gray (null-) plane is f(x, λ) = 0. For a constant value of λ, the dynamics is determined by the black curve. The fixed points are determined by the intersection with the null plane (green); the point in the front is the stable fixed point, while the further point is the unstable fixed point. When changing λ toward λc = 0, the two fixed points merge at the saddle-node bifurcation (m, λc) (green). The normal form of the saddle node is obtained by expanding f(x, λ) to the lowest order around the point (m, λc), noting that f(m, λc) = 0, \({\partial }_{x}f{(x,\lambda )}_{(x,\lambda )=(m,{\lambda }_{c})}=0\) and \({\partial }_{\lambda }f{(x,\lambda )}_{(x,\lambda )=(m,{\lambda }_{c})} < 0\) (see Fig. 8):

$$f(x,\, \lambda )\, \approx \, \frac{1}{2}{\partial }_{{xx}}f{(x,\, \lambda )}_{(x,\lambda )=(m,\, {\lambda }_{c})}\times {(x-m)}^{2}+{\partial }_{\lambda }f{(x,\, \lambda )}_{(x,\, \lambda )=(m,{\lambda }_{c})}\\ \times (\lambda -{\lambda }_{c})=-A{(x-m)}^{2}-\tilde{\lambda },$$

(12)

where \(A=-\frac{1}{2}{\partial }_{{xx}}f{(x,\,\lambda )}_{(x,\lambda )=(m,{\lambda }_{c})}\) and \(\tilde{\lambda }=-{\partial }_{\lambda }f{(x,\lambda )}_{(x,\lambda )=(m,{\lambda }_{c})}\times (\lambda -{\lambda }_{c})\). This is the normal form for the saddle-node bifurcation. Thus, close to the bifurcation point, the stable steady state is

$${x}_{0}=m+\sqrt{-\tilde{\lambda }/A}.$$

(13)

In order to see that this is indeed the case for the AMOC transition also in comprehensive climate models, Fig. 2 is adapted from the model intercomparison study34. The steady state curves obtained are from simulations, with a very slowly changing control parameter (freshwater forcing). The top panel shows ocean-only models, while the bottom panel shows atmosphere-ocean models. The curves are, even away from the transition, surprisingly well fitted by Eq. (13). Note that for some models, the transition happens before the critical point, as should be expected from noise-induced transitions. Note also that the data has been smoothed such that increasing variance close to the transition is not visible. This observation strongly supports the assumption of a saddle-node bifurcation, while it also shows that (m, λc) (black dots) are quite different between models, thus calling for reliable determination from observations.