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.

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 [4], are connected to the bimodal nature of the AMOC [5,6].Numerous climate model studies show a hysteresis behaviour, 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 inter-model spread is large and the critical threshold is poorly constrained [10,11,12].
When complex systems 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 [13].Beside a decline of the AMOC before the critical transition, there are EWSs, statistical quantities, which also change before the tipping happens.These are critical slow down (increased auto-correlation) and, from the Fluctuation-Dissipation Theorem, increased variance in the signal [14,15,16].The latter is also termed "loss of resilience", especially in the context of ecological collapse [17].The two EWSs are statistical equilibrium concepts.Thus, using them as actual predictors of a forthcoming transition, rely 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 [18].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 [19,20,21].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.Secondly, there is an ongoing positive linear trend in the SST related to global warming, which is also not related to the circulation.This is compensated for by subtracting 2× the global mean (GM) SST anomaly (small seasonal cycle removed).This differs slightly from ref. [11], where 1× the GM SST was subtracted.The factor 2 is the optimal value for the polar amplification [22] obtained by calibrating to recent direct measurements [23] (supplementary text S6).Fig. 1b shows the SG and the GM SSTs obtained from the Hadley Centre Sea Ice and Sea Surface Temperature data set (HadISST) [24].Fig. 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,25] that this and similar AMOC indices show significant trends in the mean, the variance and the autocorrelation, indicating earlywarning of a shutdown of the AMOC.However, a trend in the EWSs within a limited period of observation could be a random fluctuation within a 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.We find that the transition where a control parameter reaches the critical value, is most likely to occur around 2057 with 95% confidence interval 2034-2128.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 [26], thus we do not assume the control parameter known.Boers [25] assumes the global mean temperature T to represent the control parameter.T increases roughly linear with time since ∼ 1920 (Fig. 1d).All we assume here is that the AMOC is in an equilibrium state prior to a change towards the transition.The simplest uninformed assumption is that the change is sufficiently slow and that the control parameter approach the (unknown) critical value linearly with time.This assumption is confirmed a posteriory by a close fit to the observed AMOC fingerprint.

Modeling and detecting the critical transition
Denote the observed AMOC fingerprint by x(t) (Fig. 1e).We model it by a stochastic process X t which, depending on a control parameter λ < 0, is in 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 X t , but can assume an effective dynamics, which, with λ sufficiently close to the critical value λ c = 0, can be described by the stochastic differential equation (SDE): where µ = m + |λ|/A is the stable fix point of the drift, A is a time scale parameter, B t is a Brownian motion and σ 2 scales the variance.Disregarding the noise, this is the normal form of the co-dimension one saddle-node bifurcation [13] (supplementary text S5).The square-root dependence of the stable state: µ − m ∼ √ λ c − λ 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 t 0 , λ(t) begins to change linearly towards λ c = λ(t c ) = 0:  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. [27].
where Θ[t] is the Heaviside function and τ r = t c − t 0 > 0 is the ramping time up to time t c , where the transition eventually will occur.Time t c is denoted the tipping time, however, it can happen earlier due to a noise-induced tipping.As the transition is approached, the risk of a noise-induced tipping (n-tipping) prior to t 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, λ) = 1 − exp(−t/τ n (λ)), with mean waiting time τ n (λ) = (π/ |λ|) exp(8|λ| 3 2 /3σ 2 ) (supplementary text S4).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, the shorter is the window T w required for detecting a difference.Conversely, close to the transition critical slow down decreases the number of independent points within a window, thus calling for a larger window for a 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 [28].A Taylor expansion around the mean µ(λ) yields the approximation where µ(λ) = m + |λ|/A and α(λ) = 2 |λ|/A is the inverse correlation time.For fixed λ the process is stationary, with mean µ, variance γ 2 = σ 2 /2α and one-lag autocorrelation ρ = exp(−α∆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: α = − log ρ/∆t, σ 2 = 2αγ 2 , A = α/2(µ − m) and λ = (σ 2 /4γ 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 [29]).The uncertainty is expressed through the variances of the estimators γ2 and ρ obtained from the observations within a time window T w .Before the ramping where the process is stationary, the uncertainties can be made arbitrarily small by observing over a long time window.We may therefore assume that ρ 0 = exp(−α 0 ∆t) and γ 2 0 = σ 2 /2α 0 are known, where α 0 = 2 |λ 0 |/A and λ 0 is the baseline value before t 0 .Detection of an EWS at some chosen confidence level q (such as 95% or 99%) requires one of the estimates γ2 or ρ for a given window to be statistically different from the baseline values, which depends 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 in the steady state.The period τ 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 T w (q, α) to detect a change from baseline in EWSs at the given confidence level q is given in closed form in the next section, (eq.( 7) for variance and eq.( 8) for autocorrelation).The involved time scales are summarized in Fig. 3a, where the required window size T 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 (blue curve).With T w = 50yrs, 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 75yrs 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 = 110yrs, indicated by the time labels on top of Fig. 3a.Ten of these realizations are shown in Fig. 3b on top of the stable and unstable branches of fixed points of model (1) (the bifurcation diagram).Fig. 3c (d) shows the variance (autocorrelation) calculated from the realizations within a running 50yrs window (shown in Fig. 3c).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, The time remaining before t c is shown on top of the plot.The red and orange curves shows the time window, T w , needed in order to detect 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 significant increase in autocorrelation requires a longer data window than detecting a significant increase in variance.With T w = 50yrs (red dot-dashed line) an increase in variance can only be detected at the 95% confidence level after the red curve is below the 50yrs level.The blue curve shows the mean waiting time for a noise-induced transition, when this becomes shorter than the 50yrs level the EWS is no longer relevant, due to n-tipping occurring before t 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).Panel b shows ten model realizations of the ramped approach to t c , notice a few n-tippings prior to t c .The black (black dashed) curve is the stable (unstable) fixed point of the model.Panel c shows the increased variance as EWS: Black line is the pre-ramping steady state value, while dashed lines are the two-sigma uncertainty range for calculating variance within the 50yr 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 Panel b, while the green band corresponds to the green band in Panel a.The thin blue lines are the same obtained from simulating 1000 realizations.Panel d is the same as Panel c but for the autocorrelation, where now the green band is narrower, corresponding to T win (ac) being smaller than the window size.
respectively.The corresponding light blue curves are obtained numerically from the 1000 realizations.The green band in Fig. 3c corresponds to the green band in Fig. 3a and shows where early warning is possible in this case.

Statistics of Early Warning Signals
The asymptotic variances of the estimators are Var(μ) where T w = n∆t is the observation window.
The question is then how large T w needs to be to detect a statistically significant increase compared to the baseline values γ 2 0 and ρ 0 .For a given estimate γ2 , the estimated difference from the baseline variance is and the estimated difference from the baseline autocorrelation is Since the two EWSs, γ 2 and ρ, are treated on an equal footing, in the following we let ψ denote either of the estimators ( 13) or ( 14), the standard error is s( ψ) = Var( ψ) 1/2 (eq.( 4)) and ∆ denotes either of the two estimated differences (5) or (6).The null hypothesis is that λ = λ 0 , or equivalently α = α 0 .The null distribution of ψ 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 ∆ < qs( ψ) 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 ∆(t) > q(s( ψ(t)) + s(ψ 0 )), which, solved for T w gives for variance: and for autocorrelation, Substituting α(t) = 2 A|λ(t)|, provides the time window T w needed to detect an EWS at time t with large probability.

Predicting a forthcoming collapse of the AMOC
The AMOC fingerprint shown in Fig. 1e (replotted in Fig. 4a) shows an increased variance, γ 2 , and autocorrelation, ρ, plotted in Fig. 4b and c as functions of the mid-point of a 50yrs running window, i.e., the EWS obtained in 2020 is assigned to year 1995.The estimates leave the confidence band of the baseline values (pink area) around year 1970.This is not the estimate of t 0 , 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 approaching 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.The first method is moment-based and uses the variance and autocorrelation estimates within the running windows.2. The second method uses 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 t 0 .In general, MLE is statistically the preferred method of choice giving the most accurate results with the lowest estimation variance.1. Moment estimator of the tipping time Within the running window, we obtain the parameters α(t) (Fig. 4d) and σ 2 (Fig. 4e) of the linearized dynamics, eq. ( 3).Then we obtain Aλ(t) from σ 2 and γ 2 (t) (Fig. 4f), using that Aλ(t) = (σ 2 /4γ 2 (t)) 2 .This is consistent with a linear ramping of λ(t) beginning from a constant level λ 0 at a time t 0 .By sweeping t 0 from 1910 to 1950 and T w from 45 to 65 yrs, we obtain Aλ 0 and τ r from least square error fit to the data.This shows a single minimum at t 0 = 1924 and T w = 55yrs (Fig. 5b).Setting t 0 = 1924, we obtain t c from a linear fit (regressing λ on t) from the crossing of the x-axis (λ c = 0).This is shown in Fig. 4f (red line).This yields −Aλ 0 = 2.34 year −2 and τ r = 133 years.Thus, the tipping time is estimated to be in year 2057, shown in Fig. 4f.Since we have only obtained the combined quantity Aλ = (σ 2 /4γ 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 The estimates are shown by the red curves in Fig. 4a-f.The red dot in Fig. 4a is the tipping point and the dashed line in Fig. 4b is the asymptote for the variance.
With the parameter values completely determined, the confidence levels are calculated: The two-sigma levels around the baseline values of the EWS are shown by purple bands in Fig. 4b and c.Thus, both EWSs show increases beyond the two-sigma level from 1970 and onwards.

Maximum likelihood estimator of the tipping time We use approximate MLE on 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 t 0 the 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 t 0 are used.After time t 0 , 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 [30] (supplementary text S2).
The optimal fit is t 0 = 1927 and t c = 2069 with a 95% confidence interval 2034-2128 obtained by bootstrap (see below).These estimates are close to the estimates obtained by the moment method.
Uncertainty in the estimate of the tipping time 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 repeat 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.
From the tipping times estimated on each simulated data set, the probability density function (PDF) (Fig. 4f, yellow histogram) is obtained.The median is t c = 2066 and the 95% confidence interval is 2034 − 2128.The small discrepancy in median 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 (supplementary text S3) were calculated for the data.These are plotted in Fig. 5f 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.

Summary
We have provided a novel robust statistical analysis to quantify the uncertainty in observed EWSs for a forthcoming critical transition.The confidence depends on how rapid 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, a Kendall's τ test.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 2057.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 steadily change of the control parameter towards the collapse of the AMOC (i.e.reduce temperature increase and fresh water input through ice melting into the North Atlantic region).As a collapse of the AMOC has strong societal implications [31], it is important to monitor the flow and EWS from direct measurements [32,33,34].

Supplementary text S1 Maximum likelihood estimators of the Ornstein-Uhlenbeck process
To obtain eq. ( 4) we need the maximum likelihood estimator (MLE) of the approximate model.The approximate model is an Ornstein-Uhlenbeck (OU) process, defined as the solution to the equation This is a Gaussian process with well-known properties [29,35].The variance is γ 2 = σ 2 /2α and the ∆t-lag autocorrelation is ρ = e −α∆t .The likelihood function of the parameters given observations (x 0 , x 1 , . . ., x n ) is the product of the transition densities where θ = (µ, ρ, γ 2 ).Here, see [29,35] for details.The likelihood is the joint probability of the observed data viewed as a function of the parameters of the statistical model, in this case discrete observations from the Ornstein-Uhlenbeck process.Considering the observed sample as fixed, the likelihood is a function of the parameters.
The likelihood principle states that all the information about the parameter θ is given in the likelihood function.The maximum likelihood estimator is the value of θ which maximizes the probability of observing the given sample.In practice, the maximum of the likelihood function is found by taking the derivative with respect to the parameters (the score) and equate it to zero (the likelihood equation).For further details about likelihood theory, see any textbook in mathematical statistics.The maximum likelihood estimators (MLEs) derived from eqs. ( 10) and ( 11) are the symbol ˆindicates an estimator.These are obtained as follows.The score function is the vector of derivatives of the log-likelihood function with respect to the parameters.The MLE is given as solution to the likelihood equations ∂ θ k log L n (θ) = 0, where θ k is either µ, ρ or γ 2 .The score function is whose zeros provide the MLEs in equations ( 12)-( 14).It requires that The Fisher Information I of the MLEs equals minus the expectation of the Hessian H of the log-likelihood function.For the OU log-likelihood, the elements of H are given by where C i , i = 1, . . ., 5, are deterministic constants that will disappear when taking expectations.Using that , we obtain the Fisher Information The inverse of the Fisher Information provides the asymptotic covariance matrix, The diagonal elements provide the asymptotic variances of µ, ρ and γ 2 , respectively.

S2 Estimator of the tipping time
The process is given as solution to and we wish to estimate the parameters θ = (A, m, λ 0 , τ r , σ) from observations (x 0 , x 1 , . . ., x n ) before time t 0 and observations (y 0 , y 1 , . . ., y n ) after time t 0 , of process X t defined by (15).This equation cannot be explicitly solved, and the exact distribution is not explicitly known.A standard way to solve this is approximating 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 [30].Instead we use a two-step procedure: First we estimate α 0 = 2 A|λ 0 |, µ 0 = m + |λ 0 |/A and σ 2 from the stationary part before time t 0 , using estimators ( 12) -( 14), where α 0 = − log(ρ)/∆t and σ 2 = 2α 0 γ 2 .This yields estimates λ 0 (A) = −α 2 0 /4A and m(A) = µ 0 − α 0 /2A as a function of parameter A and the estimated parameters.The two remaining parameters A and τ r are then estimated from the data after time t 0 , where we no longer can use the OU process, since the linear approximation breaks down when the tipping point is approached.Simplifying by assuming that λ is constant between observations, i.e., piecewise constant and jumping every month where new AMOC observations are available, we obtain transition densities that are non-linear transformations of Gaussian densities, making the inference problem tractable as follows.We use a pseudo-likelihood induced by the Strang splitting scheme, shown to be robust for highly non-linear models [30].Consider the two subsystems where α(λ) = 2 A|λ| and µ(λ) = m + |λ|/A.The drift of subsystem (17) is the Taylor expansion of the drift in eq. ( 15) to first order around the fixed point µ(λ) and is an OU process, of which we know the distribution and the likelihood (see S1). Eq. ( 18) is a deterministic equation with the non-linear part, which solution is also known.We obtain the following two flows: where ).The Strang splitting [30] then approximates by which is defined for all x > µ(λ t ) − 2/A∆.Since we are only interested in simulating the process up to the time where X t crosses the separatrix between the two attractors, which happens for x < m, we require that m > m + |λ t |/A − 2/A∆ ≥ m + |λ 0 |/A − 2/A∆, i.e., ∆ < 2/ A|λ 0 | = 4/α 0 .This is always fulfilled, since ∆ = 1/12 and α 0 is estimated to be less than 4.
The transition density ( 19) is a nonlinear transformation of a Gaussian random variable, leading to the pseudo-loglikelihood function (up to a constant) where ∆/2 (y i−1 ) + µ(λ ti−1 )(1 − e −α(λt i−1 )∆ ), see [30] for details.The first two terms in (20) are the standard terms from a Gaussian likelihood, the last term originates from the non-linear transformation.Estimates of parameters A, τ r are then obtained by minimizing (20).Since division by A enters the calculations of λ 0 and m and thus the pseudo-likelihood, estimates are sensitive to small values of A. We therefore regularize the optimization problem by adding a penalization term on small values of A. The term −pn(1/A − 1) is added to (20) for A < 1, where p ≥ 0 is a penalization parameter determined by cross-validation on simulated data sets by minimizing the mean squared distance between the estimated ramping time on each data set to the value of the ramping time used in the simulation.The optimal value was p = 0.004.
The parameter estimates are found numerically by minimizing − log L x (θ).For this we apply the optimizer optim in R, using the Nelder-Mead algorithm.
Confidence intervals are obtained by parametric bootstrap: 1000 repetitions of the model are simulated with the estimated parameters, and on each synthetic data set, parameters are estimated.The empirical quantiles of the 1000 estimates thus obtained are used to construct confidence intervals.

S3 Model control
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 observation x i follows some distribution function F i, θ for the estimated parameter values θ.If this is true, then u i = F i, θ (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.

S4 Noise induced tipping
The drift term in eq. ( 1) is the negative gradient of a potential, f ( For λ < 0, the drift has two fixed points, m ± |λ|/A.The point m + |λ|/A is a local minimum of the potential V (x, λ) and is stable, whereas m − |λ|/A is a local maximum and unstable.The system thus has two basins of attraction separated by m − |λ|/A, with a drift towards either m + |λ|/A or −∞ dependent on whether X t > m − |λ|/A or X t < m − |λ|/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 X t will be drawn towards −∞.
Due to the noise, the process can escape into the tipped state by crossing over the potential barrier ∆ . Assume X t to be close to m + |λ|/A at some time t, i.e., in the normal state.The escape time will asymptotically (for σ → 0) follow an exponential distribution such that 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 (λ) is [36,37]: 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.

S5 Normal form of the saddle-node bifurcation
Consider the general dynamical equation where x is a variable and λ is a (fixed) parameter.A point x 0 with f (x 0 , λ) = 0 is a fix point or steady state.A fix point is stable/unstable if ∂ x f (x, λ) x=x0 is negative/positive, thus the fix point is attracting/repelling.If f (x, λ) is not a linear function of x, multiple steady states may exist.A saddle-node bifurcation occur when changing the control parameter λ through a critical value λ c a stable and an unstable fix point merge and disappears.The situation is shown in the figure, where the blue surface is f (x, λ), while the grey (null-) plane is f (x, λ) = 0.For a constant value of λ the dynamics is determined by the black curve.The fix points are determined by the intersection with the null-plane (green), the point in the front is the stable fix point, while the further point is the unstable fix point.When changing λ towards λ c = 0, the two fix points merge at the saddle-node bifurcation (m, λ c ) (green).The normal form of the saddle-node is obtained by expanding f (x, λ) to lowest order around the point (m, λ c ), noting that f (m, λ c ) = 0, ∂ x f (x, λ) (x,λ)=(m,λc) = 0 and ∂ λ f (x, λ) (x,λ)=(m,λc) < 0 (see Fig 6): where A = − 1 2 ∂ 2 x 2 f (x, λ) (x,λ)=(m,λc) and λ = −∂ λ f (x, λ) (x,λ)=(m,λc) × (λ − λ c ).This is the normal form for the saddle-node bifurcation.Thus, close to the bifurcation point the stable steady state is 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 [27].The steady state curves obtained are from simulations, with a very slowly changing control parameter (freshwater forcing).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. ( 25).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 calls for reliable determination from observations.The function f (x, λ) near a saddle point where a stable and an unstable fix point merge at a saddle-node bifurcation.For λ < λ c there are two fix points (green) where the black curve pass through the null-plane, f = 0 (grey).The point in front is the stable fix point, while the point in the back is the unstable fix point.The red curve of fix points is the bifurcation curves, with the stable branch in front and unstable branch in the back.The purple curve is f (x, λ c ) which touch the null-plane in one point (x 0 , λ c ).At this point it is seen that ∂ x f = 0 and ∂ λ f < 0, indicated by the dark green tangents to the surface.

S6 The AMOC proxy
The Ceasar et al. proxy is the mean SST over the SG region subtracted the global mean in order to compensate for global warming on top of the change in the AMOC.The "translation" from the proxy SST temperature and AMOC flow is 0.26 SV/K (Caesar et al. (2018), 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 Frajka-Williams et al. (2019).This shows a drop of 3 SV in that period.Minimizing the difference between the proxy SST SG -A SST GM and this more direct measurement with respect to A we get A = 1.95 ≈ 2 rather than A= 1 used by Caesar et al.The orginal and our calibrated proxies are shown in Fig. 7.

Figure 1 :
Figure 1: Panel a shows the Subpolar gyre (SG) region (black contour) on top of the HasISST SST reconstruction for Dec. 2020.The SG region SST has been identified as an AMOC fingerprint [19].Panel b shows full monthly record of the SG SST together with the global mean (GM) SST.Panels c and d show the SG and GM anomalies, which are the records subtracted the monthly mean over the full record.Panel e shows the AMOC fingerprint proxy, which is here defined as the SG anomaly minus twice the GM anomaly, compensating for the polar amplified global warming.

Figure 2 :
Figure 2: The steady state curves from climate model simulations of the North Atlantic Deep Water (NADW), with a very slowly changing control parameter (freshwater forcing).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.[27].

Figure 3 :
Figure3: Panel a shows time scales involved in the critical transition ramping the control parameter λ from λ 0 = −2.82 to λ c = 0, with a ramping time τ r = 110yrs and σ 2 = 0.29.These parameters are obtained as best estimates from the HasISST data.The time remaining before t c is shown on top of the plot.The red and orange curves shows the time window, T w , needed in order to detect 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 significant increase in autocorrelation requires a longer data window than detecting a significant increase in variance.With T w = 50yrs (red dot-dashed line) an increase in variance can only be detected at the 95% confidence level after the red curve is below the 50yrs level.The blue curve shows the mean waiting time for a noise-induced transition, when this becomes shorter than the 50yrs level the EWS is no longer relevant, due to n-tipping occurring before t 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).Panel b shows ten model realizations of the ramped approach to t c , notice a few n-tippings prior to t c .The black (black dashed) curve is the stable (unstable) fixed point of the model.Panel c shows the increased variance as EWS: Black line is the pre-ramping steady state value, while dashed lines are the two-sigma uncertainty range for calculating variance within the 50yr 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 Panel b, while the green band corresponds to the green band in Panel a.The thin blue lines are the same obtained from simulating 1000 realizations.Panel d is the same as Panel c but for the autocorrelation, where now the green band is narrower, corresponding to T win (ac) being smaller than the window size.

Figure 4 :
Figure 4: Panel a shows the SST anomaly (identical to Figure 1e) together with best estimate model of the steady state approaching a critical transition.Panels b and c show variance and autocorrelation calculated within running 50yr windows, similar to Figure 3c and d.The two-sigma level (indicated by the purple band) is obtained using the model to estimate the time varying α (Panel d) and σ 2 (Panel e) from the data.Panel f shows the best estimate for t c .The yellow histogram is the probability density for t c obtained by maximum likelihood estimates (see Methods).

Figure 5 :
Figure 5: With parameters obtained from the data, a set of 1000 realizations of the model are used in a bootstrap study to assess the uncertainty on parameters.Panels a-d are probability densities for t c , λ 0 , m and −A × m.Red crosses are the values obtained from the AMOC fingerprint data.The 95% confidence intervals are indicated by orange lines.The critical time t c is 2057, and the 95 % confidence interval is 2034-2128.Panel e shows the mean square error in fitting the ramping as a function of window size T w and time of initiating ramping, t 0 .A unique minimum is found for T w = 55 yrs and t 0 = 1924.Panel d shows the QQ-plot of residuals from the model, if points fall close to a straight line (black line) the model fits the data well.

Figure 6 :
Figure6: The function f (x, λ) near a saddle point where a stable and an unstable fix point merge at a saddle-node bifurcation.For λ < λ c there are two fix points (green) where the black curve pass through the null-plane, f = 0 (grey).The point in front is the stable fix point, while the point in the back is the unstable fix point.The red curve of fix points is the bifurcation curves, with the stable branch in front and unstable branch in the back.The purple curve is f (x, λ c ) which touch the null-plane in one point (x 0 , λ c ).At this point it is seen that ∂ x f = 0 and ∂ λ f < 0, indicated by the dark green tangents to the surface.

Figure 7 :
Figure 7: In the SST AMOC proxy the compensation for global warming and polar amplification is done by subtracting the global SST from the SG SST.By calibrating by the MOC z AMOC proxy (red curves) the optimal AMOC proxy is SST SG -2 SST GM .