Assessing AMOC stability using a Bayesian nested time-dependent autoregressive model

The Atlantic Meridional Overturning Circulation (AMOC) is a major climate element subject to possible ongoing loss of stability. Recent studies have found evidence of a gradual weakening in circulation, including early warning signals (EWSs), such as increased fluctuations and correlation time of the system, which are both known to be indicators of a possible forthcoming tipping point. To assess these changes in statistical behavior, we propose a robust and general statistical model based on a second-order autoregressive process with time-dependent parameters. This allows for the statistical changes from increased external variability and destabilization to be accounted for separately. We estimate the time evolution of the correlation parameters using a hierarchical Bayesian modeling framework, which also yields uncertainty quantification through the posterior distribution. To assess possible changes in AMOC stability, we apply the model to an AMOC fingerprint proxy based on the subpolar gyre and the global mean temperature anomaly. We find statistically significant EWSs, which suggests that AMOC is indeed undergoing a loss of stability and is getting closer to a tipping point. The methodology developed in this study is made publicly available as an extension of the R-package INLA.ews.

1 Introduction

The Atlantic Meridional Overturning Circulation (AMOC) is a key driver of Earth's climate, responsible for the transport of heat and salt across the Atlantic Ocean (Rahmstorf, 1995). As part of the global thermohaline circulation, the AMOC plays a central role in maintaining the current climate equilibrium. It is widely believed that the AMOC is a multi-stable system, capable of existing in multiple stable modes, most notably a strong mode, which is currently dominant, and a weak or collapsed mode (Stommel, 1961; Lenton et al., 2008). This nonlinear behavior implies that the AMOC may undergo abrupt transitions between states when critical thresholds are crossed. Paleoclimate evidence supports the idea that abrupt shifts in AMOC strength have contributed to major climate events during the last glacial period, such as the Dansgaard–Oeschger events (Vettoretti and Peltier, 2016; Boers et al., 2018).

These dynamics have led to the identification of the AMOC as a potential “tipping element” in the Earth system, i.e., a subsystem that could undergo a critical transition due to anthropogenic forcing (Lenton et al., 2008). Climate models suggest that continued greenhouse gas emissions and the resulting increase in freshwater input from Greenland Ice Sheet melt, precipitation and river discharge could push the AMOC toward such a tipping point (Wood et al., 2019; Hawkins et al., 2011; Weijer et al., 2019). This behavior exhibits hysteresis, meaning that once a tipping threshold is passed, the AMOC may not return to its original state even if the perturbation is reversed.

Recent observational and modeling studies have intensified concerns. Although early models suggested a low probability of collapse within the 21st century (Masson-Delmotte et al., 2021), more recent simulations reveal a wider range of possible responses, raising concerns that risks might be underestimated (Gong et al., 2022). This discrepancy is partly due to model biases, notably in representing freshwater forcing and feedback (Liu et al., 2017). Evidence is also emerging from real-world observations. Studies have documented a significant weakening trend in the AMOC over the 20th century (Caesar et al., 2018), and recent statistical analyses have detected early warning signals of reduced stability (Boers, 2021; Ditlevsen and Ditlevsen, 2023). These findings suggest that the AMOC may be approaching a critical threshold.

A weakening of the AMOC would have profound and potentially irreversible consequences, including disrupting weather patterns, altering precipitation systems and potentially triggering cascading effects on other climate components (Stouffer et al., 2006; Jackson et al., 2015; Lenton et al., 2008). In light of this, there is an urgent need to monitor the resilience of the system and improve our understanding of the processes that drive its potential loss of stability. To anticipate such changes, studies have focused on using critical slowing down theory, stating that when a system is approaching a tipping point, its recovery from small perturbations becomes progressively weaker. This phenomenon, called an early warning signal (EWS), can be characterized by increased variance and autocorrelation, which can be used as statistical indicators of approaching critical transitions.

To detect these statistical changes, a common approach is to start from the linear approximation of a dynamical system with white noise around some stable fixed point x_s, giving

$$\begin{array}{}\text{(1)}& \mathrm{d}x\left(t\right)=-\mathit{\lambda }\left(x\left(t\right)-x_{\mathrm{s}}\right)\mathrm{d}t+\mathit{\sigma }dB\left(t\right),\end{array}$$

where λ is the restoring rate and dB(t) is a white noise process. This linearization is recognized as the Langevin stochastic differential equation, which has the following solution:

$$\begin{array}{}\text{(2)}& x\left(t\right)=x_{\mathrm{0}}+\underset{-\mathrm{\infty }}{\overset{t}{\int }}g(t-s)dB\left(s\right),\end{array}$$

where g(t−s) is a Green's function defined by

$$\begin{array}{}\text{(3)}& g\left(t\right)=\left\{\begin{array}{ll}\exp (-\mathit{\lambda }t),& x\ge \mathrm{0}\\ \mathrm{0},& x<\mathrm{0}\end{array}\right..\end{array}$$

This form of x(t) is also referred to as an Ornstein–Uhlenbeck (OU) process. When discretized, this process yields a first-order autoregressive (AR) process:

$$\begin{array}{}\text{(4)}& x_t=\mathit{\varphi }{x}_{t-\mathrm{1}}+{\mathit{\epsilon }}_t,{\mathit{\epsilon }}_t\sim \mathcal{N}\left(\mathrm{0},{\displaystyle \frac{\mathrm{1}-{\mathit{\varphi }}^{\mathrm{2}}}{\mathrm{2}\mathit{\lambda }}}{\mathit{\sigma }}^{\mathrm{2}}\right),\end{array}$$

with variance $\mathrm{Var}\left(x_t\right)={\mathit{\sigma }}^{\mathrm{2}}/\left(\mathrm{2}\mathit{\lambda }\right)$ and lag-one autocorrelation parameter $\mathit{\varphi }=\exp (-\mathit{\lambda }\mathrm{\Delta }t)$.

With this model, EWSs are detected through an increase in the autocorrelation or variance. However, Boers (2021) showed that these indicators can be biased if the system is driven by external noise that itself has increasing autocorrelation or variance, leading to false positive alarms. To account for such bias, Boettner and Boers (2022) and Morr and Boers (2024) suggest that the OU process of Eq. (2) should be driven by correlated noise rather than white noise. After discretization, the resulting process yields an AR(1) process that is driven by another AR(1) process. Hence, the discretization is similar to Eq. (4), except that the white noise process ε_t is replaced by an AR(1) process:

$$\begin{array}{}\text{(5)}& v_{t+\mathrm{1}}=\mathit{\rho }{v}_t+{\mathit{\sigma }}_v{\mathit{\xi }}_t,\end{array}$$

where ρ represents the correlation parameter of the noise, σ_v is a scaling parameter and

$$\begin{array}{}\text{(6)}& {\mathit{\xi }}_t\sim \mathcal{N}\left(\mathrm{0},{\displaystyle \frac{\mathrm{1}-{\mathit{\varphi }}^{\mathrm{2}}}{\mathrm{2}\mathit{\lambda }}}\right)\end{array}$$

is a white noise process. This model encompasses the original AR(1) model in Eq. (4) when ρ=0 and, as showcased in Boers (2021), comprehends cases in which external noise is also correlated, preventing bias in the estimation of the parameter ϕ. Consequently, an increasing ϕ will act as a more reliable indicator for detecting EWSs, as it will no longer be affected by rising external variation.

Climate systems that are prone to tipping, such as the AMOC, are often driven by some external forcing. For the AMOC, the freshwater forcing from Greenland melts acts like a bifurcation parameter, as freshwater inputs can disturb the salinity and the temperature of the AMOC, potentially pushing the system closer to its tipping point (Wood et al., 2019). To incorporate forcing into our model, we use a similar approach as in Myrvoll-Nilsen et al. (2024) and Myrvoll-Nilsen et al. (2020), where the dynamical system is represented by

$$\begin{array}{}\text{(7)}& \mathrm{d}x\left(t\right)=-\mathit{\lambda }x\left(t\right)+F\left(t\right)\mathrm{d}t+U\left(t\right)\mathrm{d}t,\end{array}$$

where F(t) represents the forcing and U(t), as before, represents an OU process. The solution of this equation can be expressed as the sum of one forced component and one noise component:

$$\begin{array}{}\text{(8)}& x\left(t\right)=\mathit{\nu }\left(t\right)+\mathit{\xi }\left(t\right).\end{array}$$

Here, the noise component, ξ(t), is represented by a nested AR(1) process described previously, and the forced component, ν(t), is expressed by

$$\begin{array}{}\text{(9)}& \mathit{\nu }\left(t\right)={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\lambda }\left(t\right){\mathit{\kappa }}_f}}}\underset{\mathrm{0}}{\overset{t}{\int }}F\left(s\right)e^{-\mathit{\lambda }\left(t\right)(t-s)}\mathrm{d}s,\end{array}$$

with κ_f being a scaling parameter. This model allows EWSs to be detected while accounting for the influence of external forcing on the system's dynamics.

Most studies detect EWSs using sliding windows to obtain estimates of the variance and correlation for each window. This approach requires selecting an appropriate window length, which introduces a fundamental compromise. A shorter window provides a more accurate representation of the system's momentary state, but the limited number of data points can reduce the reliability of the statistical estimates. In contrast, a longer window improves the robustness of these estimates by incorporating more data, but it does so at the cost of responsiveness, as it averages information over a broader timescale and may fail to capture short-term fluctuations effectively. Determining the optimal window length is thus a critical but challenging task, as it should ideally balance estimation accuracy with the ability to reflect rapid changes in the system's evolution. Myrvoll-Nilsen et al. (2024) propose an alternative model-based approach that eliminates the need for this choice. Instead of relying on a fixed window length, the correlation parameter is assumed to evolve over time according to a predefined linear structure. This assumption enables a hierarchical Bayesian model formulation, enabling the use of well-established computational techniques to infer the parameters of the linear structure. Furthermore, Myrvoll-Nilsen et al. (2024) adopts a Bayesian framework that offers the additional advantage of providing uncertainty quantification in the form of posterior distributions, making the analysis more robust and interpretable.

In this paper, we build upon the hierarchical Bayesian framework developed by Myrvoll-Nilsen et al. (2024) to integrate the nested AR(1) model proposed by Morr and Boers (2024) and Boettner and Boers (2022). This extension helps mitigate false alarms caused by correlated noise and eliminates the need for sliding time windows, while benefiting from the advantages of a Bayesian modeling framework. This approach is then applied to an AMOC fingerprint in order to assess its potential loss of stability.

The paper is organized as follows. Section 2 outlines our methodology for the Bayesian modeling framework, including details on how inference can be obtained efficiently. In Sect. 3, we evaluate our model's accuracy and reliability on simulated data, assess the robustness to false alarms under increasing external variability, and benchmark its performance on real data against existing approaches. In Sect. 4, we use our Bayesian framework to identify EWSs in an AMOC fingerprint, using different detrending strategies. Further discussion and conclusions are provided in Sect. 5.

2 Bayesian modeling

We assume that the observations, $\mathit{y}=(y_{\mathrm{1}},\mathrm{\dots },y_n{)}^{\top }$, are expressed by

$$\begin{array}{}\text{(10)}& \mathit{y}=\mathit{\mu }+\mathit{x},\end{array}$$

where the forcing response, $\mathit{\mu }=({\mathit{\mu }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\mu }}_n{)}^{\top }$, is expressed by

$$\begin{array}{}\text{(11)}& {\mathit{\mu }}_t={\mathit{\sigma }}_f\left(t\right)\sum _{s=\mathrm{0}}^{t}F\left(s\right)e^{-\mathit{\lambda }\left(t\right)(t-s+\mathrm{0.5})}\mathrm{d}s\end{array}$$

and the correlated time-dependent noise, $\mathit{x}=(x_{\mathrm{1}},\mathrm{\dots },x_n{)}^{\top }$, is given by a nested AR(1) process

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& x_{t+\mathrm{1}}=\mathit{\varphi }{x}_t+v_{t+\mathrm{1}}\\ & v_{t+\mathrm{1}}=\mathit{\rho }{v}_t+{\mathit{\sigma }}_v{\mathit{\xi }}_t.\end{array}\end{array}$$

To model the evolution of the autocorrelation parameters, we assume that they both change linearly in time, i.e.,

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& \mathit{\varphi }\left(t\right)=a_{\mathit{\varphi }}+b_{\mathit{\varphi }}t,\mathrm{0}\le t\le \mathrm{1},\\ & \mathit{\rho }\left(t\right)=a_{\mathit{\rho }}+b_{\mathit{\rho }}t,\mathrm{0}\le t\le \mathrm{1}.\end{array}\end{array}$$

These are expressed by unknown parameters a_ϕ, b_ϕ, a_ρ and b_ρ, which are estimated by fitting the model to observed data. Early warning signals due to critical slowing down are characterized through the evolution of ϕ(t), while potential changes in external variability are captured by the latent component $\mathit{v}=(v_{\mathrm{1}},\mathrm{\dots },v_n{)}^{\top }$. Separating these signals prevents false alarms, as discussed by Boers (2021).

To obtain robust uncertainty estimates, we adopt a Bayesian framework for parameter estimation, similar to Myrvoll-Nilsen et al. (2024). Given the hierarchical nature of the model, where y is modeled in terms of μ and x, which are themselves governed by hyperparameters $\mathit{\theta }=(a_{\mathit{\varphi }},b_{\mathit{\varphi }},a_{\mathit{\rho }},b_{\mathit{\rho }},{\mathit{\sigma }}_v,{\mathit{\sigma }}_f)$, a latent Gaussian model formulation provides a natural and efficient framework for Bayesian inference. Both components of the model, μ and x, depend on the parameters a_ϕ and b_ϕ through $\mathit{\lambda }\left(t\right)=-\log \mathit{\varphi }\left(t\right)$. This dependency introduces a challenge for obtaining reliable inference, as the parameters may be difficult to estimate independently. We therefore choose to model the sum $\mathit{\eta }=\mathit{\mu }+\mathit{x}$ as a single component. The latent Gaussian model formulation is defined in three stages as follows.

1. The first stage defines the likelihood of the model, which is assumed to be conditionally independent given the latent components. Because the observations y are captured by the latent component $\mathit{\eta }=({\mathit{\eta }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\eta }}_n{)}^{\top }$, we model y as a Gaussian distribution with mean η and negligible observation noise, σ_y≈0, effectively setting y≈η, i.e.,

   $$\begin{array}{}\text{(14)}& \begin{array}{rl}\mathit{\pi }(\mathit{y}\mid \mathit{\eta },\mathit{\theta })& =\prod _{k=\mathrm{1}}^n\mathit{\pi }\left(y_k\mid {\mathit{\eta }}_k,\mathit{\theta }\right)\\ & =\prod _{k=\mathrm{1}}^n{\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}}}\exp \left(-{\displaystyle \frac{{\left(y_k-{\mathit{\eta }}_k\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_y^{\mathrm{2}}}}\right).\end{array}\end{array}$$

2. The second stage defines the prior distribution for the latent field η, given the parameter θ. This component is assigned a multivariate Gaussian prior distribution with mean vector μ and a covariance matrix corresponding to the nested AR(1) process above with time-dependent ϕ(t) and ρ(t), i.e.,

   $$\begin{array}{}\text{(15)}& \mathit{\pi }(\mathit{\eta }\mid \mathit{\theta })={\mathcal{N}}_n(\mathit{\mu },\mathbf{\Sigma }).\end{array}$$

   Because η follows a nested AR(1) process, which is equivalent to an AR(2) process (Morr and Boers, 2024), its precision matrix, $\mathit{Q}={\mathbf{\Sigma }}^{-\mathrm{1}}$, is a sparse matrix of bandwidth 2. This property enables the use of computationally efficient algorithms that substantially reduce the overall computational cost.

3. The final stage defines the prior distributions for the model parameters:

   $$\begin{array}{}\text{(16)}& \mathit{\pi }\left(\mathit{\theta }\right)=\mathit{\pi }\left(b_{\mathit{\varphi }}\right)\mathit{\pi }(a_{\mathit{\varphi }}\mid b_{\mathit{\varphi }})\mathit{\pi }\left(b_{\mathit{\rho }}\right)\mathit{\pi }\left(a_{\mathit{\rho }}\mid b_{\mathit{\rho }}\right)\mathit{\pi }\left({\mathit{\sigma }}_v\right)\mathit{\pi }\left({\mathit{\sigma }}_f\right).\end{array}$$

   We assign uniform prior distributions on b_ϕ, (a_ϕ∣b_ϕ), b_ρ and (a_ρ∣b_ρ) and gamma distributions on ${\mathit{\kappa }}_v=\mathrm{1}/{\mathit{\sigma }}_v^{\mathrm{2}}$ and ${\mathit{\kappa }}_f=\mathrm{1}/{\mathit{\sigma }}_f^{\mathrm{2}}$. Note that because we assume that both $\mathrm{0}<\mathit{\varphi }\left(t\right)<\mathrm{1}$ and $\mathrm{0}<\mathit{\rho }\left(t\right)<\mathrm{1}$, then the parameter space of a_ϕ and a_ρ depends on the current state of b_ϕ and b_ρ, respectively.

The joint posterior distribution for the parameters is given by

$$\begin{array}{}\text{(17)}& \mathit{\pi }(\mathit{x},\mathit{v},\mathit{\theta }\mid \mathit{y})={\displaystyle \frac{\mathit{\pi }(\mathit{y}\mid \mathit{x},\mathit{v},\mathit{\theta })\mathit{\pi }(\mathit{x},\mathit{v}\mid \mathit{\theta })\mathit{\pi }\left(\mathit{\theta }\right)}{\mathit{\pi }\left(\mathit{y}\right)}},\end{array}$$

where π(y) is the marginal likelihood, or evidence, of y. In particular, we are interested in the marginal posterior distribution of b_ϕ, which can be obtained by integrating out the other parameters, ${\mathit{\theta }}_{-b_{\mathit{\varphi }}}$, and latent variables

$$\begin{array}{}\text{(18)}& \mathit{\pi }\left(b_{\mathit{\varphi }}\mid \mathit{y}\right)=\int \mathit{\pi }(\mathit{\theta },\mathit{x},\mathit{v}\mid \mathit{y})\mathrm{d}{\mathit{\theta }}_{-b_{\mathit{\varphi }}}\mathrm{d}\mathit{x}\mathrm{d}\mathit{v}.\end{array}$$

Because solving this integral analytically is often impossible to do in practice, the common approach is to instead approximate it using sampling-based approaches like Markov chain Monte Carlo (MCMC) methods (Robert et al., 1999). However, because the precision matrix of the latent Gaussian field is sparse, we can employ a number of computationally efficient algorithms for fast Bayesian inference. Specifically, we evaluate all marginal posterior distributions using the framework of integrated nested Laplace approximations (INLA) (Rue et al., 2009, 2017), which is particularly suited for these types of models. INLA is available as an R package at https://www.r-inla.org/ (last access: 11 April 2025) and presents a computationally superior alternative to MCMC. Because our model requires specific implementation using the custom modeling framework of R-INLA, we have decided to make the code available as a new feature in the user-friendly R-package INLA.ews, originally developed for the model described in Myrvoll-Nilsen et al. (2024). The nested time-dependent AR(1) model can be fitted by prompting

$$\begin{array}{rl}& \mathtt{\text{results}}<-\mathtt{\text{inla.ews}}(\mathtt{\text{data=y}},\mathtt{\text{forcing=z}},\\ & \mathtt{\text{model="ar2"}}).\end{array}$$

A more extensive demonstration of the package can be found in Myrvoll-Nilsen et al. (2024).

3 Assessing model accuracy and robustness

To evaluate the accuracy and robustness of the proposed time-dependent nested AR(1) model, we perform three tests. Two tests use simulated data, one from the nested AR(1) model and one from stochastic differential equations representing dynamical systems with and without loss of stability. These tests are both based on 500 independent simulated time series of length 150, matching the length of the AMOC fingerprint time series used in Sect. 4. Finally, we fit our model to a real data example – the Dansgaard–Oeschger events – in order to compare our model's results with existing methodologies. All tests are made using R-INLA with the prior distributions described in the previous section.

3.1 Model accuracy on simulated data

For the first test, we assess whether the model can recover known parameter values when fitted to simulated data generated from the same nested time-dependent AR(1) process. For each simulation, the slope parameters b_ϕ and b_ρ are independently drawn from a uniform distribution $\mathcal{U}(-\mathrm{0.9},\mathrm{0.9})$. Thereafter, the intercepts a_ϕ and a_ρ are drawn from uniform distributions with boundaries that depend on the simulated slope parameters, ensuring that the resulting ϕ(t) and ρ(t) remain within the interval (0, 1) for all time steps.

Figure 1 Results of the accuracy test for n_r=500 simulated time series of length n=150. Panels (a) and (b) show the posterior marginal mean estimated by INLA for ϕ and ρ, respectively. The blue line shows the true b used in the simulation. Panels (c) and (d) show the estimated posterior probability of the slope being positive against true values of ϕ and ρ, respectively. The vertical red line separates the true positive and negative values, while the horizontal one indicates the probability threshold 0.95 used here to determine statistical significance.

We compute the root mean square error (RMSE) between the true slope values and their marginal posterior means, $\widehat{b_{\mathit{\varphi }}}$ and $\widehat{b_{\mathit{\rho }}}$. We find the RMSE to be 0.145 for b_ϕ and 0.278 for b_ρ. We then assess whether the model reliably infers the sign of the slopes by comparing the marginal posterior probabilities $P(b_{\mathit{\varphi }}>\mathrm{0}\mid \mathit{y})$ and $P(b_{\mathit{\rho }}>\mathrm{0}\mid \mathit{y})$ to the true value of the slopes. We consider the slope for ϕ(t) and ρ(t) to be significantly positive if the posterior probabilities exceed the threshold $\mathrm{1}-\mathit{\alpha }=\mathrm{0.95}$, i.e., $P(b_{\mathit{\varphi }}>\mathrm{0}\mid \mathit{y})>\mathrm{0.95}$ and $P(b_{\mathit{\rho }}>\mathrm{0}\mid \mathit{y})>\mathrm{0.95}$, respectively. If an estimated $\widehat{b_{\mathit{\varphi }}}$ is classified as positive, given the $P(b_{\mathit{\varphi }}>\mathrm{0}\mid \mathit{y})>\mathrm{0.95}$ threshold, we count it as a true positive if the true slope is also positive, i.e., b_ϕ>0. On the other hand, if the true slope is negative, we count the estimate as a false positive. Conversely, if $P(b_{\mathit{\varphi }}>\mathrm{0}\mid \mathit{y})\le \mathrm{0.95}$, we count it as a true negative if the true slope is also negative and as a false negative if b_ϕ>0. We also count the classifications based on the estimated $\widehat{b_{\mathit{\rho }}}$, but these are of secondary interest. The sensitivity and specificity are computed by

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& \mathrm{Sensitivity}={\displaystyle \frac{\mathit{\#}\mathrm{True}\mathrm{Positives}}{\mathit{\#}\mathrm{True}\mathrm{Positives}+\mathit{\#}\mathrm{False}\mathrm{Negatives}}},\\ & \mathrm{Specificity}={\displaystyle \frac{\mathit{\#}\mathrm{True}\mathrm{Negatives}}{\mathit{\#}\mathrm{True}\mathrm{Negatives}+\mathit{\#}\mathrm{False}\mathrm{Positives}}}.\end{array}\end{array}$$

For b_ϕ, the model achieves a sensitivity of 87.7 % and a specificity of 99.8 %. For b_ρ, the sensitivity is 72.2 %, and the specificity is 99.4 %. The results from this test are summarized in Table 1 and illustrated in Fig. 1. Repeating the test with different prior distributions similar to Myrvoll-Nilsen et al. (2024) did not show significant changes, suggesting that the model is robust to the choice of prior distributions.

Table 1 Summary of statistics from Fig. 1 (top) and Fig. 2 (bottom). Results of accuracy tests are for simulated time-dependent nested AR(1) processes showing the root mean square error (RMSE) of the estimates of b_ϕ and b_ρ given true values of the simulations (blue lines in Fig. 1a and b). We also show the sensitivity and specificity (expressed as percentages) for both parameters (bottom). Results of robustness tests are for simulated tipping and non-tipping processes. We show here the RMSE of the estimates of b_ρ given true simulated values (blue lines in panels b and d). Sensitivity and specificity are presented as percentages for each process.

Figure 2 Results of the robustness test for n_r=500 simulated time series of length n≈150. Panels (a) and (b) show the posterior marginal mean estimated by INLA from the non-tipping simulations for b_ϕ and b_ρ, respectively, plotted against the true value of b_ρ from the correlated noise. The blue line in (b) shows the true b_ρ used in the simulation. Panels (c) and (d) are similar plots for the tipping simulations. In panels (a) and (c), blue dots are associated with a statistical significance for the EWS indicator b_ϕ to be positive, i.e., $P(b_{\mathit{\varphi }}>\mathrm{0}\mid \mathit{y})>\mathrm{0.95}$, while red dots mean no statistical significance, i.e., $P(b_{\mathit{\varphi }}>\mathrm{0}\mid \mathit{y})<\mathrm{0.95}$.

3.2 Robustness to false alarms under autocorrelated external variability

In the second test, we evaluate the ability of the model to reliably distinguish genuine early warning signals from changes driven solely by correlated external variability. To do so, we simulate data from two stochastic differential equations. The first represents a system approaching a tipping point, and the second remains stable but is influenced by time-dependent autocorrelated noise. This setup follows the example in Boers (2021). The tipping process is expressed by

$$\begin{array}{}\text{(20)}& \dot{x}\left(t\right)=-x^{\mathrm{3}}+x-T+v\left(t\right),\end{array}$$

where T increases linearly from −1 to 1 and v(t) is a time-dependent AR(1) process with parameters drawn in the same way as in the first test. The non-tipping process is generated by

$$\begin{array}{}\text{(21)}& \dot{x}\left(t\right)=-\mathrm{5}x+v\left(t\right),\end{array}$$

with the same structure for v(t) as in Eq. (20). Each simulation is run until the tipping point is reached (for the tipping system) or for 150 time units (for the non-tipping system), resulting in time series of approximately 150 points. The same inference methodology and classification thresholds are used here, with the distinction that an early warning signal is said to be detected when $P(b_{\mathit{\varphi }}>\mathrm{0}\mid \mathit{y})>\mathrm{0.95}$. For the tipping processes, the model correctly detected an EWS signal in 471 out of 500 simulations, corresponding to a sensitivity of 94.2 %. For the non-tipping processes, 23 out of 500 simulations were incorrectly classified as an EWS, resulting in a specificity of 95.4 %. These results, presented in Table 1 and Fig. 2, indicate that the model effectively identifies true loss of stability while maintaining a low false positive ratio, even in the presence of strongly autocorrelated noise.

To illustrate the benefits of accounting for the bias introduced by correlated noise, we also test the time-dependent AR(1) model proposed by Myrvoll-Nilsen et al. (2024), which does not separate external noise autocorrelation from loss of stability. We apply this model to the same set of simulated non-tipping processes. As illustrated in Fig. 3, this model yields posterior marginal mean estimates of b_ϕ that correlate with the values of b_ρ, rather than remaining centered around zero, as expected in the absence of a true loss of stability. In contrast, our nested AR(1) model maintains stable estimates of b_ϕ across all values of b_ρ, as shown in Fig. 2a, demonstrating its robustness to external noise.

Figure 3 Results of the robustness test from the Myrvoll-Nilsen et al. (2024) model applied on the same data as in Fig. 3a. The posterior marginal mean is plotted against the true value of b_ρ from the correlated noise. The blue line is a linear regression on the data, showing the drift of the estimates of b_ϕ. Turquoise dots are associated with a statistical significance for the EWS indicator b_ϕ to be positive, while red dots mean no statistical significance.

The simpler AR(1) model also exhibits a significantly higher rate of false positives, misclassifying 116 out of 500 simulations as tipping events, an approximately 400 % increase compared to the nested AR(1) model. Moreover, the false positive rate increases systematically with higher values of b_ρ, further highlighting the susceptibility of this model to bias from autocorrelated noise. In contrast, false detections in the nested AR(1) model are evenly distributed across all simulations. These results emphasize the importance of explicitly modeling the correlated noise structure when assessing stability in time series data.

Overall, these tests demonstrate that the proposed methodology reliably recovers the evolution of autocorrelation parameters, performs well in detecting EWSs, and is robust to prior assumptions and to structured stochastic external variability not linked to loss of stability.

3.3 Benchmarking on real data: Dansgaard–Oeschger events

Finally, we evaluate the performance of our model on real-world data associated with well-studied critical transitions. By comparing results, we can assess how much our model agrees or disagrees with existing approaches. Specifically, we use our model to analyze abrupt Greenland warmings known as Dansgaard–Oeschger (DO) events (Dansgaard et al., 1993; Johnsen et al., 1992), which represent rapid climate fluctuations that occurred during the last glacial period, where the temperature over Greenland and the North Atlantic region increased by up to 16.5 °C within a few decades (Kindler et al., 2014). DO events are often considered the archetypal example of tipping point crossings in the climate. As such, they present a natural benchmark for evaluating different EWS approaches.

For our analysis, we pair the δ¹⁸O proxy data from the Northern Greenland Ice Core Project (NGRIP) (North Greenland Ice Core Project members, 2004; Gkinis et al., 2014; Ruth et al., 2003) with the corresponding age provided by the Greenland Ice Core Project 2005 (GICC05) (Vinther et al., 2006; Rasmussen et al., 2006; Andersen et al., 2006; Svensson et al., 2008). The data are available at https://www.iceandclimate.nbi.ku.dk/data (last accessed: 5 August 2025). The model is fitted to segments preceding the 17 most recent DO events. The selected segments are highlighted in Fig. 4.

Figure 4 The 17 most recent Dansgaard–Oeschger (DO) events (vertical black lines) in the NGRIP δ¹⁸O record plotted against the GICC05 chronology. Early warning signals are estimated by fitting the model to the Greenland stadial periods (black segments) of the data preceding each DO event.

Table 2 Table comparing the posterior probability of a positive slope $P(b_{\mathit{\varphi }}>\mathrm{0}\mid \mathit{y})$ from fitting the nested AR(1) model to the different Dansgaard–Oeschger events using a second-order polynomial detrending approach. These results are compared with the probability of a positive slope $P(b>\mathrm{0}\mid \mathit{y})$ found by Myrvoll-Nilsen et al. (2024) and p-values obtained from Boers (2018) and Rypdal (2016). Events with a statistically significant positive slope are highlighted in bold.

Figure 5 AMOC fingerprint proxy from 1870 to 2020, similar to Ditlevsen and Ditlevsen (2023), using the yearly averaged subpolar gyre sea-surface temperature anomaly minus twice the global mean anomaly obtained from the Hadley Centre Sea Ice and Sea Surface Temperature (HadISST) dataset (Rayner et al., 2003).

Whether or not DO events are induced solely by noise or if they are indeed approaching a bifurcation point is currently debated (Ditlevsen et al., 2007; Hummel et al., 2024). There is therefore no ground truth as to which, if any, DO event should exhibit an EWS. However, several studies report detection of an EWS before some of the first 17 DO events (Rypdal, 2016; Boers, 2018). We compare the results of our model with these studies using a setup similar to Myrvoll-Nilsen et al. (2024) by using a second-order polynomial detrending of the data and considering $P(b_{\mathit{\varphi }}>\mathrm{0}\mid \mathit{y})>\mathrm{0.95}$ as a detection of EWSs. This comparison is illustrated in Table 2 and shows that our model suggests, similarly to Myrvoll-Nilsen et al. (2024), that some specific event shows signs of critical slowing down in line with the results of Boers (2018) and Rypdal (2016). Specifically, Table 2 shows that our results corroborate the five EWSs found by Myrvoll-Nilsen et al. (2024) while identifying one more EWS for the 13th event. Moreover, these results corroborate the EWS found for the 11th event by Boers (2018) and the fifth and ninth events by Rypdal (2016); our results also show EWSs for the second and 13th events, similarly to these two studies.

4 Detecting early warning signals in an AMOC fingerprint

We now apply the time-dependent nested AR(1) model to an AMOC fingerprint similar to the one used by Ditlevsen and Ditlevsen (2023) shown in Fig. 5. This fingerprint is constructed as the sea-surface temperature (SST) anomaly in the subpolar gyre region, averaged annually, minus twice the global mean SST anomaly to compensate for the polar amplification effects under global warming. Several studies have suggested that this proxy is a suitable indicator of AMOC strength (Caesar et al., 2018; Jackson and Wood, 2020; Latif et al., 2019), especially because direct observations are available only from 2004 onward. The use of such a proxy is therefore necessary to examine longer-term trends and detect potential early warning signals.

Figure 6 Integrated Central-West Greenland runoff from 1871 to 2013.

Table 3 Summary of statistics from Figs. 7 and 8 showing posterior marginal means of $\widehat{b_{\mathit{\varphi }}}$, probability of $\widehat{b_{\mathit{\varphi }}}$ being positive, posterior marginal means of $\widehat{b_{\mathit{\rho }}}$ and marginal log-likelihood for the three models used here. Results from the models introduced in Myrvoll-Nilsen et al. (2024) are also shown for comparison purposes. Events with a statistically significant positive slope are highlighted in bold.

As the fingerprint exhibits significant drift, it must first be detrended to satisfy the zero-mean assumption of the model. In principle, this trend could be extracted using knowledge of the system's underlying physical processes, but such information may be unavailable, incomplete or inaccurate. To address this, we consider two different detrending strategies. In the first, we rely solely on statistical assumptions and remove the trend using either a linear or second-order polynomial fit. In the second approach, we incorporate physical information by including an explanatory variable in the model, following the structure described in Eq. (9). Specifically, we use the integrated Central-West Greenland (iCWG) surface melt shown in Fig. 6 as a covariate. The iCWG represents the cumulative surface melt across several years, based on the CWG melt stack from Trusel et al. (2018), and is used to capture the influence of freshwater forcing on AMOC stability.

For each model, we compare the posterior marginal mean estimate of the slope parameter b_ϕ along with the posterior probability that the slope is positive. Model fit is assessed using the marginal log-likelihood. The full set of results is presented in Table 3. The fitted trends and time evolutions of ϕ(t) for the linear and polynomial detrending approaches are shown in Fig. 7, while the estimated response function to the iCWG forcing and associated ϕ(t) evolution are shown in Fig. 8. Among the different model configurations, the version incorporating iCWG forcing provides the best fit to the data as measured by the model likelihood. In all three detrending strategies, the model identifies statistically significant EWSs. These results provide further evidence for the presence of EWSs for the AMOC, consistent with the findings of Boers (2021), who also found statistically significant EWSs using a slightly different nested AR(1) process with a window-based estimation methodology applied to a similar proxy for AMOC strength; the global mean temperature is only subtracted once in their study. Our results also corroborate those found by Ditlevsen and Ditlevsen (2023), who reported similar EWSs using the same proxy but applied an AR(1) model with a window-based approach.

Figure 7 Panels (a) and (b) show the AMOC fingerprint (black) with the posterior marginal mean (blue) and 95 % credible intervals (red) of the fitted trends. Panels (c) and (d) show the evolution in time of the correlation parameter ϕ(t) (blue) used as an indicator of EWSs and the 95 % credible intervals (red) with an estimated probability of a positive slope $P(b_{\mathit{\varphi }}>\mathrm{0}\mid \mathit{y})$.

Figure 8 Panel (a) shows the AMOC fingerprint (black) from 1870 to 2013 to match the time span of the forcing data with the posterior marginal mean (blue) and 95 % credible intervals (red) of the estimated system's response function to forcing. Panel (b) is a plot of the evolution in time of the correlation parameter ϕ(t) (blue) and 95 % credible intervals (red) with an estimated probability of a positive slope $P(b_{\mathit{\varphi }}>\mathrm{0}\mid \mathit{y})$.

5 Conclusions

This study investigates the stability of the Atlantic Meridional Overturning Circulation (AMOC) by proposing a time-dependent extension of the nested autoregressive AR(1) model introduced by Morr and Boers (2024) and Boers (2021). The primary objective of this model is to enhance the reliability of early warning signals (EWSs) by minimizing false positives. This is achieved through the decomposition of the observed signal into two distinct components: ρ(t), which captures time-dependent external variability, and ϕ(t), which reflects changes in the internal dynamics associated with system stability. By isolating these effects, the model aims to identify more accurately early signs of destabilization. Following the approach of Myrvoll-Nilsen et al. (2024), we assume a linear temporal dependence for both ρ(t) and ϕ(t), estimating their respective slope parameters within a hierarchical Bayesian framework. This statistical approach allows us to incorporate prior information and quantify the uncertainty of the EWSs through the posterior distributions of the parameters. The performance of the model is first evaluated using both simulated and real data, demonstrating both high estimation accuracy and robustness against false detections of ongoing destabilization.

The methodology is applied to a proxy for the AMOC fingerprint. In order to meet stationarity assumptions, we consider various detrending techniques, including linear and second-order polynomial detrending, and incorporate a forcing component based on the integrated meltwater runoff from Central-West Greenland. Across all model configurations, we find statistically significant early warning signals. This is consistent with prior findings in the literature and supports the hypothesis of a possible ongoing destabilization of the AMOC.

While assuming a linear structure for ϕ(t) has proven effective for detecting EWSs, we emphasize that the model proposed here should not be interpreted as a comprehensive or mechanistic representation of the underlying physical processes governing the AMOC. Despite its success in identifying early signs of destabilization, the model is limited in its ability to forecast the future trajectory of the system or predict the timing of a potential tipping point. Addressing these limitations would require a more flexible modeling approach, potentially involving a nonlinear or nonparametric structure for the correlation parameters, which lies beyond the scope of the present work.

Although our analysis has focused on a specific proxy of the AMOC fingerprint, the proposed methodology is generalizable and can be adapted to study the stability of other critical climate components, such as the Greenland Ice Sheet, Arctic sea ice or the Amazon rainforest. To facilitate wider use and reproducibility, we have extended the existing R package INLA.ews to incorporate our methodological advancements. This software provides a user-friendly interface for implementing our approach, leveraging the computational efficiency of the INLA framework for Bayesian inference.