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 system^1,2,3. In recent years, model studies and paleoclimatic reconstructions indicate that the strongest abrupt climate fluctuations, the Dansgaard-Oeschger events⁴, are connected to the bimodal nature of the AMOC^5,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 bifurcations^7,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 constrained^10,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)¹². 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 collapse¹³. There are, however, model biases toward overestimated stability of the AMOC, both from tuning to the historic climate record¹⁴, poor representation of the deep water formation¹⁵, salinity and glacial runoff¹⁶.

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 dynamics¹⁷. 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 signal^18,19,20. The latter is also termed “loss of resilience”, especially in the context of ecological collapse²¹. 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 measurements²². 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 AMOC^23,24,25.

a Subpolar gyre (SG) region (black contour) on top of the Hadley Centre Sea Ice and Sea Surface Temperature data set (HadISST) SST reconstruction for December 2020. The SG region SST has been identified as an AMOC fingerprint²³.b Full monthly record of the SG SST together with the global mean (GM) SST.c,d SG and GM anomalies, which are the records subtracted the monthly mean over the full record.e AMOC fingerprint proxy, which is here defined as the SG anomaly minus twice the GM anomaly, compensating for the polar amplified global warming.

Full size image

Figure 1b shows the SG and the GM SSTs obtained from the Hadley Centre Sea Ice and Sea Surface Temperature data set (HadISST)²⁶. 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 reported^11,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 constrained²⁸; thus, we do not assume the control parameter known. Boers²⁷ assumes the global mean temperatureT to represent the control parameter. AlthoughT 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 CO₂ concentration, does, in fact, increase close to linearly with time in the industrial period²⁹. 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 byx(t) (Fig. 1e). We model it by a stochastic processX_t, 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λ = λ₀. We are uninformed about the dynamics governing the evolution ofX_t but can assume effective dynamics, which, withλ sufficiently close to the critical valueλ_c = 0, can be described by the stochastic differential equation (SDE):

where\(m=\mu -\sqrt{|\lambda|/A}\) andμ is the stable fixed point of the drift,A is a time scale parameter,B_t is a Brownian motion andσ² scales the variance. Disregarding the noise, this is the normal form of the co-dimension one saddle-node bifurcation¹⁷ (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.

Top panel shows ocean-only models, while bottom panel shows atmosphere-ocean models. The curves are, even away from the transition, surprisingly well fitted by Eq. (1) (black thin curves). The bifurcation points are indicated with black circles. Note that for some models, the transition happens before the critical point, as should be expected from noise-induced transitions. The colored circles show the present-day conditions for the different models. Adapted from Rahmstorf et al.³⁴.

Full size image

At timet₀,λ(t) begins to change linearly towardλ_c = λ(t_c) = 0:

where Θ[t] is the Heaviside function andτ_r = t_c − t₀ > 0 is the ramping time up to timet_c, where the transition eventually will occur. Timet_c is denoted the tipping time; however, an actual tipping can happen earlier than t_c due to a noise-induced tipping. As the transition is approached, the risk of noise-induced tipping (n-tipping) prior tot_c 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,T_w, 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 windowT_w 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 − T_w/2, t + T_w/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 process³⁰. A Taylor expansion around the fixed pointμ(λ) yields the approximation

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α 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μ,γ² andρ the parameters of Eq. (1) are determined:\(\alpha=-\log \rho /{{\Delta }}t\),σ² = 2αγ²,A = α/2(μ − m) and\(\lambda={({\sigma }^{2}/4{\gamma }^{2})}^{2}/A\). Closed form estimators forμ, γ² andρ are obtained from the observed time series within a running window by maximum likelihood estimation (MLE) (Supplementary text S1, see also ref.³¹).

The uncertainty is expressed through the variances of the estimators\({\hat{\gamma }}^{2}\) and\(\hat{\rho }\) obtained from the observations within a time windowT_w. 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 levelq (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,t_ac, 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α(λ) andq-percentile, the required time windowT_w(q, α) to detect a change from baseline in EWSs at the given confidence levelq 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 sizeT_w 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,t_noise, (blue curve). WithT_w = 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λ₀ toλ_c, i.e., the ramping timeτ_r. A set of 1000 realizations has been simulated withλ₀ = −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λ = λ₀, 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.

a Bifurcation diagram, the blue curve is the stable state (mean),μ, and the dashed red curve is the unstable branch representing the edge state,x_e, both as functions of the control parameterλ. Atλ_c the system has a saddle-node bifurcation.b Linear ramping (Eq. (2)), where the system resides in a stationary state withλ = λ₀ until timet₀, from whereλ increases linearly with time and reach the critical valueλ_c at timet_c.c Effective dynamics (Eq. (1)) forλ = λ₀ indicated by the potentialU(x) = A(x−m)³/3 + λx (yellow) along the same color dashed line in (a). The red parabola is the quadratic potential corresponding to the approximating linear Ornstein–Uhlenbeck (OU) process (Eq. (3)). The curve in gray is a realization of the OU process. The vertical short red bar shows the autocorrelation time,t_ac, for the process. The long vertical red bar symbolizes the time window,T_w for calculating the early-warning signals (EWS).d As (c) forλ₁ close toλ_c (purple dashed line in (a)). It can be seen directly that the gray realization has an increased autocorrelation and increased variance compared to the situation in (c).

Full size image

a Time scales involved in the critical transition ramping the control parameterλ fromλ₀ = −2.7 toλ_c = 0, with a ramping timeτ_r = 133 years andσ² = 0.3. These parameters are obtained as best estimates from the Hadley Centre Sea Ice and Sea Surface Temperature data set (HadISST). The time remaining beforet_c is shown on top of the plot. The red and orange curves show the time window,T_w, needed in order to detect the increase in variance (red) and autocorrelation (orange) above the pre-ramping values at the 95% confidence level. Close to the bifurcation point, the (quasi-)stationarity approximation becomes less valid, which is indicated by the dashed part of the two curves. It is seen that detecting a significant increase in autocorrelation requires a longer data window than detecting a significant increase in variance. WithT_w = 50 years (red dot-dashed line), an increase in variance can only be detected at the 95% confidence level after the red curve is below the 50-year level. The blue curve shows the mean waiting time for a noise-induced transition; when this becomes shorter than the 50-year level, the early-warning signal (EWS) is no longer relevant due to n-tipping occurring beforet_c. Thus the range of time where an EWS can be applied is indicated by the green band (limited by the crossings of the red and blue curves with the size of the window).b Ten model realizations of the ramped approach tot_c, notice a few n-tippings prior tot_c. The black (black dashed) curve is the stable (unstable) fixed point of the model.c Increased variance as EWS: Black line is the pre-ramping steady state value, while the dashed lines are the two-standard error uncertainty range for calculating variance within the 50-year data window. The blue and dashed blue curves are the same but for the model approaching the transition. The brown curves correspond to the ten realizations in (b), while the green band corresponds to the green band in (a). The thin blue lines are the same obtained from simulating 1000 realizations.d Same as (c) but for the autocorrelation, where now the green band is narrower, corresponding toT_win(ac) being smaller than the window size.

Full size image

Statistics of early-warning signals

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

whereT_w = nΔt is the observation window.

The question is then how largeT_w 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

and the estimated difference from the baseline autocorrelation is

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λ = λ₀, or equivalentlyα = α₀. The null distribution of\(\hat{\psi }\) is assumed to be Gaussian (confirmed by simulations). A quantileq from the standard Gaussian distribution expresses the acceptable uncertainty in measuring the statistical quantityψ. We thus get that\(\hat{{{\Delta }}} \, < \, qs(\hat{\psi })\) at theq-confidence level (95%, 99% or similar) under the null hypothesis. To detect an EWS at theq-confidence level based on measuringψ at timet, we require that\(\hat{{{\Delta }}}(t) \, > \, q(s(\hat{\psi }(t))+s({\hat{\psi }}_{0}))\), which, solved forT_w gives for variance:

and for autocorrelation,

Substituting\({\alpha }_{0}=2\sqrt{A|{\lambda }_{0}|}\) and\(\alpha (t)=2\sqrt{A|\lambda (t)|}\) provides the time windowT_w needed to detect an EWS at timet 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,γ², 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 oft₀, 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.

a Sea surface temperature (SST) anomaly (identical to Fig. 1e) together with the best estimate model of the steady state approaching a critical transition.b,c Variance and autocorrelation calculated within running 50-year windows, similar to Fig. 4c, d. The two-standard error level (indicated by the purple band) is obtained using the model to estimate the time varyingα (d) andσ² (e) from the data.f Best estimate fort_c. The yellow histogram is the probability density fort_c obtained by maximum likelihood estimates (see Supplementary Information S1 andS2).

Full size image

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 timet₀. 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σ² (Fig. 5e) of the linearized dynamics, Eq. (3), and thus alsoγ²(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 obtainAλ(t) fromσ² andγ²(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λ₀ at a timet₀. By sweepingt₀ from 1910 to 1950 andT_w from 45 to 65 years, we obtainAλ₀ andτ_r from the least square error fit to the data. This shows a single minimum att₀ = 1924 andT_w = 55 years (Fig. 6e). Settingt₀ = 1924, we obtaint_c from a linear fit (regressingλ ont) from the crossing of thex-axis (λ_c = 0). This is shown in Fig. 5f (red line). This yields −Aλ₀ = 2.34 year⁻² 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 determineA andm 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.

With parameters obtained from the data, a set of 1000 realizations of the model are used in a bootstrap study to assess the uncertainty of parameters.a–d Probability densities fort_c,λ₀,m and −A × m. Red crosses are the values obtained from the Atlantic meridional overturning circulation (AMOC) fingerprint data. The 95% confidence intervals are indicated by orange lines. The critical timet_c is 2057, and the 95% confidence interval is 2025–2095.e Mean square error in fitting the ramping as a function of window sizeT_w and time of initiating ramping,t₀. A unique minimum is found forT_w = 55 years andt₀ = 1924.f Quantile-quantile-plot of residuals from the model, if points fall close to a straight line (black line), the model fits the data well.

Full size image

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 timet₀, approximation (3) is used, where exact MLEs are available (Supplementary text S1). This provides estimates of the parametersλ₀, m as a function of parameterA, as well as the variance parameterσ². To estimateA andτ_r, the observations after timet₀ are used. After timet₀, 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 inaccurate³² (Supplementary text S2). Usingt₀ = 1924, the optimal fit is the same as the moment method,t_c = 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 〈t_c〉 = 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.

We have provided a robust statistical analysis to quantify the uncertainty in observed EWSs for a forthcoming critical transition. The confidence depends on how rapidly the system is approaching the tipping point. With this, the significance of the observed EWSs for the AMOC has been established. This is a stronger result than just observing a significant trend in the EWS by, say, Kendall’sτ test^27,33. Here we calculate when the EWS are significantly above the natural variations. Furthermore, we have provided a method to not only determine whether a critical transition will happen but also an estimate of when it will happen. We predict with high confidence the tipping to happen as soon as mid-century (2025–2095 is a 95% confidence range). These results are under the assumption that the model is approximately correct, and we, of course, cannot rule out that other mechanisms are at play, and thus, the uncertainty is larger. However, we have reduced the analysis to have as few and sound assumptions as possible, and given the importance of the AMOC for the climate system, we ought not to ignore such clear indicators of an imminent collapse.

The hysteresis simulations gathered in the model intercomparison³⁴ are equilibrium runs, for which a prediction of a future collapse does obviously not apply. Likewise, for the simulations specified in the CMIP6 experiment. It could though be relevant to evaluate our method on state-of-the-art climate model simulations with linearly ramped external forcing and different ramping speeds in order to obtain the model-specific confidence in early prediction of the collapse judged solely from the EWSs.

Though we have established firm statistical methods to evaluate the confidence in the observed EWS, we can at present not rule out the possibility that a collapse will only be partial and not lead to a full collapse of the AMOC as suggested by some models: Note in Fig. 2, the “MPM” in the top panel and the “MOM hor” in the bottom panel both seem to show only partial tipping prior to the tipping to the complete shutdown of the AMOC. This result is also found in a more recent ocean model³⁵. Furthermore, a high speed of ramping, i.e., a high speed at which the critical value of the control parameter is approached, could also increase the probability of tipping³⁶. This scenario is the case of rate-induced tipping³⁷. Even with these reservations, this is indeed a worrisome result, which should call for fast and effective measures to reduce global greenhouse gas emissions in order to avoid the steady change of the control parameter toward the collapse of the AMOC (i.e., reduce temperature increase and freshwater input through ice melting into the North Atlantic region). As a collapse of the AMOC has strong societal implications³⁸, it is important to monitor the flow and EWS from direct measurements^39,40,41.

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.²³, where 1 × the GM SST was subtracted. The “translation” from the proxy SST temperature and AMOC flow is 0.26 SV/K [ref.²³, 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 MOC_z reported in the review by ref.⁴². This shows a drop of 3 SV in that period. Minimizing the difference between the proxy SST_SG-β SST_GM and this more direct measurement with respect toβ, we getβ = 1.95 ≈ 2 rather thanβ = 1 used by ref.²³. The factor 2 is thus the optimal value for the polar amplification⁴³ obtained by calibrating to recent direct measurements⁴². The original and our calibrated proxies are shown in Fig. 7.

In the sea surface temperature (SST) AMOC fingerprint record, the compensation for global warming and polar amplification is done by subtracting the global SST (SST_GM) from the subpolar gyre (SG) SST (SST_SG). By calibrating by the MOC_z AMOC proxy (red curves), the optimal AMOC proxy is SST_SG-2 SST_GM. To ensure the robustness of our results, we have repeated the analysis by subtracting 1x (purple) and 3x SST_GM (green) and compared with the optimal 2x SST_GM subtracted (blue). The corresponding estimates for the time of the collapse are shown in the same colors: The middle vertical line is the maximum likelihood estimate of the tipping time, the box represents the 66.6% confidence interval (Intergovernmental Panel on Climate Change (IPCC) definition of “likely”), while the horizontal line represents the 99% confidence intervals.

Full size image

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 fixedt₀ = 1924.

Estimator of the tipping time and model control

The AMOC fingerprint is assumed to be observations from a processX_t given as a solution to Eqs. (1) and (2), and we wish to estimate the parametersθ = (A, m, λ₀, τ_r, σ) from observations (x₀, x₁, …, x_n) before timet₀ and observations (y₀, y₁, …, y_n) after timet₀. These equations cannot be explicitly solved, and the exact distribution ofX_t 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 models³². We estimate by a two-step procedure, approximating the stationary distribution before timet₀ 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 timet₀, using methods proposed in ref.³², see Supplementary text S2 for details.

To test the model fit, uniform residuals,u_i, i = 1, …, n were calculated for the AMOC data using the estimated parameters from the MLE method as follows. The model assumes that observationx_i 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)² + λ) withV(x, λ) = A(x−m)³/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 potentialV(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 andX_t 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}\). AssumeX_t to be close to\(m+\sqrt{|\lambda|/A}\) at some timet, i.e., in the normal state. The escape time will asymptotically (forσ → 0) follow an exponential distribution such that

whereP(t, λ) is the probability of observing an escape time shorter thant for a given value ofλ. The mean noise-induced escape timeτ_n(λ) is^44,45:

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

wherex is a variable andλ is a (fixed) parameter. A pointx₀ withf(x₀, λ) = 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. Iff(x, λ) is not a linear function ofx, 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 isf(x, λ), while the gray (null-) plane isf(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 expandingf(x, λ) to the lowest order around the point (m, λ_c), noting thatf(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):

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

The functionf(x, λ) near a saddle point where a stable and an unstable fixed point merge at a saddle-node bifurcation. Forλ < λ_c, there are two fixed points (green) where the black curve passes through the null plane,f = 0 (gray). The point in front is the stable fixed point, while the point in the back is the unstable fixed point. The red curve of fixed points is the bifurcation curve, with the stable branch in front and unstable branch in the back. The purple curve isf(x, λ_c) which touches the null plane at one point (x₀, λ_c). At this point, it is seen that ∂_xf = 0 and ∂_λf < 0, indicated by the dark green tangents to the surface.

Full size image

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 study³⁴. 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.

Data can be found in the following repository: https://doi.org/10.17894/ucph.b8f99b67-d4e6-4a2e-b518-00bddeed323b.

Computer code (Matlab and R) can be found in the following repository:https://doi.org/10.17894/ucph.b8f99b67-d4e6-4a2e-b518-00bddeed323b.

This work has received funding under the project Tipping Points in the Earth System (TiPES) from the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 820970. This is TiPES contribution #214. Furthermore, funding was provided by Novo Nordisk Foundation NNF20OC0062958; and European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 956107, “Economic Policy in Complex Environments (EPOC)”.

1. These authors contributed equally: Peter Ditlevsen, Susanne Ditlevsen.

Authors and Affiliations

1. Niels Bohr Institute, University of Copenhagen, Copenhagen, Denmark

   Peter Ditlevsen

2. Institute of Mathematical Sciences, University of Copenhagen, Copenhagen, Denmark

   Susanne Ditlevsen

Corresponding authors

Correspondence toPeter Ditlevsen orSusanne Ditlevsen.

Competing interests

The authors declare no competing interests.

Peer review information

Nature Communications thanks Chris Boulton and the other anonymous reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visithttp://creativecommons.org/licenses/by/4.0/.

Reprints and permissions