Limits to predictability of the asymptotic state of the Atlantic Meridional Overturning Circulation in a conceptual climate model

Anticipating critical transitions in the Earth system is of great societal relevance, yet there may be intrinsic limitations to their predictability. For instance, from the theory of dynamical systems possessing multiple chaotic attractors, it is known that the asymptotic state depends sensitively on the initial condition in the proximity of a fractal basin boundary. Here, we approach the problem of final-state sensitivity of the Atlantic Meridional Over-turning Circulation (AMOC) using a conceptual climate model, composed of a slow bistable ocean coupled to a fast chaotic atmosphere. First, we explore the occurrence of long chaotic transients in the monostable regime, which can mask a loss of stability near bifurcations. In the bistable regime, we explicitly construct the chaotic saddle using the edge tracking technique. Quantifying the final-state sensitivity through the maximum Lyapunov exponent and the lifetime of the saddle, we find that the system exhibits a fractal basin boundary with almost full phase space dimension, implying vanishing predictability of the second kind near the basin boundary. Our results demonstrate the usefulness of studying non-attracting chaotic sets in the context of predicting climatic tipping points, and provide guidance for the interpretation of higher-dimensional models such as general circulation models.


Introduction
Like many systems in nature, several elements of the Earth system are thought to be multistable: for a given climatic forcing, they may possess multiple competing attractors that can be reached from different initial conditions [1].Multistability has been demonstrated in rather complex physical models of the Greenland [2] and West Antarctic [3] ice sheets, the Amazon rainforest [4], the Atlantic Meridional Overturning Circulation (AMOC) [5], and even Earth as a whole [6], supported by paleoclimatic evidence for abrupt climate changes in the past [7,8].The proposed multistability of the AMOC has been intensively studied for several decades [9,10,11] and has recently attracted renewed attention due to a suggested loss of stability over the past century seen in observation-based indicators [12].The paradigm of a multistable AMOC dates back to the seminal model introduced by Stommel [13], who showed that there can be two competing ("on" and "off") states for a given freshwater forcing due to the positive salt-advection feedback.This concept has been invoked to explain qualitative changes of the ocean circulation in past climates [14,15] as well as the response of climate models to changes in the hydrological cycle [16,17,18].
Collectively, the multistable climate subsystems mentioned above are subsumed under the term tipping elements [19], as anthropogenic global warming may trigger a possibly irreversible transition (tipping) to a qualitatively different state on timescales relevant for policy [20].In turn, the concept of a "safe operating space" has been framed as keeping a control parameter (such as greenhouse gas forcing) below a critical value to avoid the transition into an undesired state of the Earth system [21].However, in addition to the classical picture of crossing a bifurcation ("bifurcation tipping"), internal variability and the rate of the forcing can also lead to noise-induced and rate-induced tipping, respectively [22], which means that in reality the safe operating space is likely more complex than what can be described with a single critical value [23,24,25].
By incorporating the dominant physical processes into low-dimensional conceptual models, dynamical systems theory has proven extremely useful in advancing our understanding of the mechanisms of tipping points in the Earth system (see [8,26,27,28,29] for recent reviews on this topic).For example, the derivation of statistical early warning indicators [30,31,32], which can be applied to observational (e.g., [12,33]) or modeled (e.g., [34]) climatic time series, has increased the prospect of anticipating that a bifurcation is approached.However, in a multistable chaotic system there are fundamental limitations to predictability of the final state even in the autonomous case.The question of whether the attractor reached from a given initial condition can be accurately predicted has been coined "predictability of the second kind" by Lorenz [35].This is in contrast to predictability "of the first kind", which refers to the ability to predict the future state of a system at a given horizon, given the knowledge of the initial conditions with finite precision.The presence of a limited predictability of the first kind is a key characteristic of chaotic systems and is usually investigated through Lyapunov analysis [36,37].
Let us consider the simpler case of bistable systems.Uncertainty on our ability to predict the final state emerges from the (unavoidably) finite observation time.Beyond a bifurcation point, one may encounter so-called "ghost attractors" -states that are not asymptotically stable but feature transient chaos with finite yet possibly very long lifetimes, which depend sensitively on the initial conditions [38].In practice, this means that the system appears to reside in a well-defined steady state until a sudden transition to the actual attractor occurs.
This phenomenon is common in many areas of physics; for example, turbulence in pipe flows is often regarded as a very long chaotic transient rather than a genuine attractor (e.g., [39]).In this case, given a finite time of observation or simulation (such as in state-of-the-art climate models, which are typically integrated for decades to centuries), one might misjudge the stability properties of the system.This is particularly relevant when a control parameter is close to the critical value separating monostable and bistable behavior.
In the parametric region where bistability is observed, a separate instance of difficulty emerges in predicting the asymptotic state from the choice of the initial condition.In a chaotic, bistable D-dimensional system, limited final state predictability is directly linked to the presence of a fractal boundary between the two basins of attraction with fractal dimension D − 1 ≤ D b < D [40]: given that an initial condition u 0 can only be determined to a precision ε, the fraction f of (a bounded region of) phase space in which the outcome is uncertain (i.e., different attractors can be reached from within u 0 ± ε) scales like where α = D − D b is the uncertainty exponent [40].In practical terms, α ≪ 1 means that decreasing the uncertainty ε only yields a very small improvement in final state predictability as given by f , and the phase space region around the boundary is essentially a "grey zone" in which the final state is almost unpredictable.The case α ≪ 1 is believed to be relatively common [41].It has been conjectured that the fractal basin boundary dimension is linked to the Lyapunov spectrum and the lifetime of the chaotic saddle [42,43], with α ≪ 1 being associated with a chaotic instability on the saddle that is fast compared to its lifetime [44].While their importance for the transient and asymptotic behavior of chaotic dynamical systems has long been recognized [38,45], both fractal basin boundaries and chaotic saddles remain understudied in the context of climatic tipping points.To our knowledge, the only study that explicitly determined the basin boundary dimension in the context of two competing climatic attractors has been by Lucarini & Bódai [46].They used a climate model of intermediate complexity to investigate the "Snowball Earth" transition [47,48], in which the ice-albedo feedback drives almost every initial condition either to a fully glaciated or an ice-free climate, and found that the basin boundary has almost full dimension, with α ≈ 0.02.Regarding the AMOC, Lohmann & Ditlevsen [25] found that in the context of rate-induced tipping (occurring before the bifurcation point is crossed [22]), the final AMOC state in a simplified ocean general circulation model (GCM) depends sensitively on the initial condition.They qualitatively linked this behavior to the presence of a fractal basin boundary but did not assess its properties explicitly.
In this study, we quantitatively explore the limits of predictability of the final AMOC state in a conceptual atmosphere-ocean model inspired by Gottwald [49] (Sec.2).The model mimics a key feature of more complex GCMs: it exhibits chaotic dynamics driven by a fast atmospheric component, while the oceanic component acts as a slow integrator.The oceanic component introduces the bistability in the system as a result of the coexistence of the AMOC "on" and "off" states.We investigate aspects related to both time and phase space that limit predictability of the second kind.First, we focus on the lifetime of chaotic transients which may effectively prevent the system from tipping on finite timescales, even when overshooting the bifurcation (Section 3).Then, building on the methods of [46], we use the edge tracking algorithm [50] to construct the chaotic saddle (also called Melancholia state by [46]) between two competing AMOC states, whose lifetime and Lyapunov exponents are directly linked to the dimension of the fractal basin boundary [42,43] (Section 4).Finally, we discuss the implications of both phenomena on the predictability of AMOC tipping and on the concept of a safe operating space (Section 5) and summarize our conclusions in Section 6.

The conceptual atmosphere-ocean model
First, we introduce the conceptual climate model which isolates the principal characteristics of atmosphere-ocean flow that we will focus on in this study -a fast, chaotic atmosphere coupled to a slow, bistable ocean component.To this end, we use a coupled configuration of the seminal conceptual models for mid-latitude atmospheric flow and the AMOC, respectively: the Lorenz-84 model [51] (L84 hereafter) and the Stommel model [13].

Model components
The L84 model for the atmospheric mid-latitude circulation is given by a set of three ordinary differential equations [51]: where x is the strength of the westerly winds (a conceptual representation of the jet stream) given by thermal wind balance, and y and z are the amplitudes of the cosine and sine phases of a superimposed traveling wave.The parameters F and G are the background values for the meridional and zonal temperature gradients, respectively, to which x and y would relax in an uncoupled setting; a and b control the internal timescales of the model and will be kept at their original values a = 0.25 and b = 4 [51].Despite its very simple representation of atmospheric flow, the model has a strong physical basis, as it can be derived by Galerkin truncation of a two-layer quasi-geostrophic model [52].It also exhibits a complicated bifurcation structure on its own [53,54], but since different states of atmospheric flow are not the focus of our study, we remain in the regime G = 1 where the L84 model only has a single chaotic attractor depicted in Fig. 2.
The Stommel model [13] is the canonical conceptual model for the bistable AMOC and describes the thermohaline flow in the North Atlantic as a function of the meridional temperature gradient T and the salinity gradient S between two boxes representing the subtropical and subpolar North Atlantic.For several decades, Stommel's model has served as a key model to interpret qualitative changes in the AMOC [16,10,11].Here, we use the version derived by Cessi [55] and follow the notation of Gottwald [49]: where ε a = t r /t d is the ratio between the restoring timescale for temperature t r and the diffusive timescale of the ocean t d , µ = t ad /t d is the ratio between the advective timescale t ad and the diffusive timescale, θ is the surface temperature gradient to which T is restored, and σ is the surface freshwater flux.The advective terms are proportional to the absolute value of the AMOC strength Ψ = T − S, where Ψ > 0 represents a thermally driven circulation like in the present-day climate, and Ψ < 0 represents a salinity-driven AMOC where the circulation would be reversed.Here, as in many previous studies, we consider the freshwater flux σ as the main control parameter that determines the AMOC regime.
The corresponding bifurcation diagram, which reveals a bistable region bounded by two saddle-node bifurcations, is shown in Fig. 1.

Coupled model equations
As in several previous studies [49,56,57] of conceptual atmosphere-ocean coupling, we couple the Stommel model to the L84 atmosphere through the temperature gradients  and the freshwater flux (Fig. 2).Here, we follow the recent formulation of coupling terms by Gottwald [49] (without the sea ice component), who introduced an explicit timescale separation ε f between the atmospheric and oceanic components, similar to earlier studies [56,57].The equations of the five-dimensional coupled model used in our study are given by: with Here, in a similar spirit as in [56,57], we use a two-way coupling: the fast atmospheric component drives the ocean via the restoring temperature and salinity/freshwater flux, while the slow ocean component couples to the atmosphere via the meridional temperature gradient F .For simplicity, we keep the zonal temperature gradient G = G 0 constant, noting that the chaotic properties of the L84 model as characterized by its maximum Lyapunov exponent already follow an intricate non-monotonic structure when only F is varied [58].
The main coupling constant is F 1 , which we set to F 1 = 0.1 following the similar model of [57], who argued that the atmosphere should be weakly coupled to the ocean.Furthermore, the timescale separation parameter between atmosphere and ocean is chosen as ε f = 3•10 −4 , which corresponds to a timescale of about 10 days for the atmosphere and 100 years for the ocean.Note that this oceanic timescale defines the scale of the dimensionless time t, such that all times given here can be read in a "unit" of centuries.Our oceanic timescale is slightly shorter than the one given by Cessi [55], but we mostly aim at a realistic order of magnitude, as a smaller ε f needs to be traded off against a higher computational cost for the edge tracking.Except for F 1 and ε f , we use the default parameters of [49] (see Appendix A for a full list of parameter definitions and values).The equations are solved numerically using a fifth-order Runge-Kutta scheme with a fixed timestep of 7.5•10 −6 , which corresponds to a timestep of 1/40 of the atmospheric timescale as used in [56].
We remark that it is possible to perform a rigorous asymptotic analysis of Eq. 4 that, in the limit ε f → 0, leads to writing the homogenized dynamics for the slow variables (T, S) in terms of a stochastic differential equation.The presence of a two-way coupling between the fast and the slow variables requires a more general approach [59] than the more common case of skew-symmetric systems, where the dynamics of the fast variables has an impact on the slow variables but not the other way around [60,61].See Appendix B for a sketch of the asymptotic multiscale analysis of the system given in Eq. 4 and a derivation of the corresponding homogenized evolution equations for the slow variables.

Chaotic transients
We first turn our attention to the limits of predictability caused by the sensitivity of the lifetime of chaotic transients to the initial condition.As an example, we select a set of initial conditions by sampling the "AMOC on" state for σ 0 = 0.926 every ∆t = 7.5.As we will show explicitly below (cf.Fig. 8), at this value of σ 0 the "AMOC on" state is globally stable.Then, we integrate trajectories from these initial conditions but with the freshwater parameter set to σ 0 = 0.932 (Fig. 3), i.e. slightly larger than for the bifurcation point L 2 of the uncoupled Stommel model, σ(L 2 ) = 0.9263.We can assume that this value is also close to the bifurcation point of the coupled model because the coupling is weak (see Fig. 8).
The initial conditions in Fig. 3 remain in the vicinity of a "ghost attractor" reminiscent of the "AMOC on" state for a long time (up to 400 time units ≈ 40 000 years) before converging to the only remaining genuine attractor, the "AMOC off" state.In this section, we systematically explore the lifetime of such chaotic transients and show that they can be observed over a wide range of σ 0 .For an ensemble of initial conditions, the lifetime of chaotic transients, denoted as τ , is expected to be exponentially distributed [62] with a mean lifetime ⟨τ ⟩, if τ is sufficiently large [63].Indeed, this is the case for a 100-member initial condition ensemble from which the trajectories in Fig. 3 were drawn, with ⟨τ ⟩ ≈ 200 time units (20 000 years).Grebogi et al. [63] then demonstrated that typically, the mean lifetime follows a power law as one approaches the critical value of the control parameter σ 0,c , such that we expect: where γ is the critical exponent.
Here, we characterize this region of transient chaos by evaluating the lifetimes of trajectories near the "ghost attractor" using an ensemble of initial conditions for different values of σ 0 .The initial conditions are sampled over t = 200 from the "off" or "on" state near the bifurcation points of the underlying Stommel model, and the results are not sensitive to the exact location of the sampling.Here and in the following, all sampling intervals are chosen to be at least several times larger than the Lyapunov timescale of the fast system (≈ 0.002 time units).Because ⟨τ ⟩ → ∞ as σ 0 → σ 0,c , we need to choose a threshold time T max after which we deem a trajectory "stable", even though it may still be a very long chaotic transient.This means that σ 0,c would be biased if determined by direct simulation, as it would depend on the choice of T max .Therefore, we look to find a range of values for σ 0 with ⟨τ ⟩ sufficiently large for the values to be exponentially distributed, but smaller than T max .
Using Eq. 5, we then determine σ 0,c and γ within this region via a weighted least-squares fit, taking into account that the parameter ⟨τ ⟩ is distributed as 2n ⟨ τ ⟩/χ 2 (2n) [64], where ⟨ τ ⟩ is equal to the sample mean of the lifetimes, n is the number of samples, and χ 2 (2n) is the chi-square distribution with 2n degrees of freedom.the lifetime and |σ 0 − σ 0,c | near both boundary crisis points.We obtain σ 0,c1 = 0.816 ± 0.008 and σ 0,c2 = 0.9273 ± 0.0009 (cf.Fig. 9b), where index "1" denotes the transition from the monostable "AMOC on" to the bistable regime, and index "2" denotes the transition from the bistable to the monostable "AMOC off" regime.These values are quite robust to the choice of the scaling law.For instance, fitting the non-power law scaling proposed by [65], we obtain σ 0,c2 = 0.9251 ± 0.0013, which agrees well with the estimate from the fit to Eq. 5.The critical exponents are γ 1 = 14 ± 3 at the first and γ 2 = 7.2 ± 1.5 at the second bifurcation.This implies that long transients can be observed over a wider range of parameter values at the first compared to the second bifurcation [45].For instance, the parameter range |σ 0 − σ 0,c | for which ⟨τ ⟩ > 100 is 0.04 at the first bifurcation and 0.005 at the second bifurcation.Nevertheless, both critical exponents are larger than those of many typical dynamical systems [45].Using the scaling law from Eq. 5 and the quantile function of the exponential distribution, we can now calculate how far the "safe operating space" of the AMOC is extended beyond the genuinely bistable regime through the existence of long transients.In other words, up to which value of σ 0 will the desirable AMOC state (in our case, the "on" state) continue to act like a stable attractor, even though the bifurcation point has been crossed?The exponential distribution underlying the transient lifetime means that this safe operating space is inherently linked to a timescale of interest and an acceptable probability of survival on the desired AMOC state (or, equivalently, a tolerable probability of tipping to the undesired state).Fig. 5 shows this probability of survival on the "on" state as a function of |σ 0 − σ 0,c | and the chosen timescale.For example, overshooting to σ 0 ≈ 0.934 would yield a 5% probability of obtaining a transient lifetime of less than 100 years (red marker in Fig. 5).On the other hand, if the probability of a "transient collapse" should be limited to 0.1% on the same timescale, the safe operating space is limited to no more than σ 0 ≈ 0.931.In this way, if the control parameter is changed sufficiently rapidly to bring the system back into the bistable regime, a transition might successfully be avoided despite "overshooting" the critical value (cf.[66]).

Chaotic saddle and fractal basin boundary
We now focus on the limits to predictability in the bistable regime.The long-lived transients make it difficult to estimate its boundaries precisely, such that we rely on the two extrapolated critical values from Sec. 3 to bound this regime (σ 0 ∈ [0.82, 0.927] in the following), within which the "AMOC on" and "AMOC off" states are expected to coexist as genuine attractors.As outlined above, in the bistable regime predictability of the asymptotic state is limited by the dimension D b of the fractal basin boundary.D b is in turn fundamentally linked to the properties of the chaotic saddle between the two attractors [42,43,44].However, the chaotic saddle is typically a set with Lebesgue measure zero that has one unstable and one stable manifold (the basin boundary) [45], and therefore cannot be obtained by direct simulation.Therefore, we start with constructing a pseudo-trajectory of the saddle (Sec.4.1) which we use to determine its lifetime (Sec.4.2) and Lyapunov exponents (Sec.4.3) before assessing the resulting basin boundary dimension (Sec.4.4).

Constructing the chaotic saddle: edge tracking algorithm
While there are different algorithms to construct a chaotic saddle approximately [38], we follow the strategy of [46] by using the edge tracking algorithm originally proposed in [67] and later re-formulated by [50].This allows us to obtain an arbitrarily long pseudo-trajectory that tracks the chaotic saddle of the conceptual atmosphere-ocean model.The edge tracking algorithm has previously been applied successfully in both high-and low-dimensional chaotic systems (e.g., [46,50,68,69,70]).The algorithm is depicted schematically in Fig. 6 and described in detail in Sects.4.2 and 5.1 of [46], such that we only summarize it briefly here.We start from two initial conditions {u 0 1 , u 0 2 } inside the two different basins of attraction.First, by applying the standard bisection method, a pair of states {u 1 and u * 2 still belong to two different basins of attraction.By construction, these two states are both within a distance δ 1 from the basin boundary.As a consequence, when the system is evolved from u i (t 0 ) ≡ u * i , they are expected to diverge along the unstable direction of the basin boundary while tracking with the flow along its stable direction.Once |u 1 (t) − u 2 (t)| ≥ δ 2 for some t > t 0 , the bisection is repeated such that the distance between the two "shadowing trajectories" u 1 (t) and u 2 (t) is again < δ 1 , and the new states are evolved further.This procedure can be repeated any number of times to obtain a pseudo-trajectory u S (t) = (u 1 (t) + u 2 (t))/2 that follows the basin boundary and eventually converges to track the chaotic saddle itself.
In the following, we set δ 2 = 0.004 and δ 1 = 0.0025, such that only one bisection step is needed per iteration.However, as in more complex models, it is not clear a priori how the norm |u 1 (t) − u 2 (t)| should be defined, as the state vector u comprises different model components with possibly different magnitudes for mean and variability.Here, we define the norm by rescaling the atmospheric variables with a scaling factor ξ < 1 (we choose ξ = 1/5000).It can be viewed as a hyperparameter of our method and does not alter the dynamics of the dynamical system.This enables us to calculate the divergence of trajectories over the full phase space (instead of just the T −S subspace), but ensures that the divergence is not dominated by the fast chaotic motion.We run the edge tracking algorithm for 2000 iterations, generating around 50 time units (≈ 5000 years) of pseudo-trajectory.
Fig. 7 shows an example result of the edge tracking algorithm for σ 0 = 0.9.While the attractors and the saddle separate clearly in both T and S as expected from the phase space structure of the underlying Stommel model, the atmospheric motion differs little between the two attractors.However, close to the saddle the atmospheric energy 1  2 (X 2 +Y 2 +Z 2 ) has longer tails, which might be due to the repelling along the unstable direction, and therefore a sign that the weak coupling does exert a certain influence on the atmosphere.Applying the edge tracking algorithm for different values of σ 0 , we can compute the unstable branch of the AMOC bifurcation diagram for the coupled atmosphere-ocean model (Fig. 8).Over most of the parameter range, our bifurcation diagram follows that of the underlying uncoupled Stommel model very closely, as expected from weak atmosphere-ocean coupling combined with a large timescale separation.Outside the bistable parameter range, we can also perform edge tracking in a region of parameter space where very long-lived chaotic transients exist; in this case, the algorithm detects an edge state lying in between the long-lived chaotic solution and the actual unique asymptotic solution (unfilled circles in Fig. 8).Note that in this regime, the approach is very similar to the one taken in turbulence [50], where one finds the unstable solution between the laminar fixed point and the long-lived turbulent state.An interesting feature of the bifurcation diagram in Fig. 8 is that the stable and unstable branches do not meet at the critical values σ 0,ci .This might put into question the extrapolated critical values from Sec. 3, but we have verified the transient nature of some of the unfilled points on the "AMOC on" branch of Fig. 8 via direct simulation (e.g., Fig. 3).In these cases, the transient lifetime is longer than the length of the edge tracking trajectory (t ≈ 50).Therefore, the discrepancy seems to arise from the particular multiscale nature of our system, which merits further investigation in future studies.

Lifetime of the saddle
We now focus on two properties of the chaotic saddle, its lifetime and its spectrum of Lyapunov exponents (LEs) [45], for different values of the freshwater forcing σ 0 in the bistable regime.Both can be obtained by sampling a large number of initial conditions from the edge tracking trajectory obtained in Sec.4.1.Since a state initialized on the trajectory u S is very close to but not precisely on the actual saddle, it will remain in a phase space region B containing the saddle for some time before converging to one of the attractors.
Defining B as the bounding box of u S after the initial spin-up, the number of remaining trajectories within B is expected to decay exponentially [45]: where N 0 is the number of sampled initial conditions (N 0 = 600 in our numerical simulations), which are all within B by construction, and ⟨τ ⟩ is the mean lifetime.In practice, this exponential scaling is not expected to hold for small and large t due to the non-uniform initial distribution in B and the finite size of the ensemble, respectively.However, we can obtain ⟨τ ⟩ via a least-squares fit of the slope of log N S (t)/N 0 against t for intermediate values of t.

Lyapunov spectrum of the saddle
Having explored the global instability of the coupled model characterized by the lifetime of the saddle, we now turn towards the local instability due to its chaotic dynamics.To this end, we compute the full spectrum of Lyapunov exponents of the saddle using the standard procedure of successive Gram-Schmidt orthonormalization [71] as implemented by [72].The Lyapunov spectrum of the saddle is approximated by choosing 2000 initial conditions on u S with the longest lifetimes in B and then averaging over the Lyapunov spectra of a subset of these individual trajectories while they remain in B. For consistency, we tested the same method for the attractors, where averaging over 2000 finite-time LEs calculated over a similar duration (t = 7.5) yielded a comparable result to letting the algorithm converge along one long trajectory.
Uncertainty estimates are 95% confidence intervals generated from bootstrap resampling of 100 trajectories each from the 2000-member ensemble.
The Lyapunov spectra {λ i }, i = 1, ..., 5 for the saddle and each of the two attractors are given in Tab. 1 for σ 0 = 0.9.For the attractors, we expect one positive LE characterizing the chaotic instability on the attracting set, a zero LE associated with the motion along the attractor, and otherwise negative LEs encapsulating the convergence of initial conditions to the attractor.Indeed, we find a large, positive maximum Lyapunov exponent (MLE) λ 1 , while λ 2 is very close to zero and all remaining LEs are negative.Similarly, the saddle possesses one positive LE (λ 1 ), one LE close to zero (λ 2 ) and three negative LEs.Note that obtaining a good numerical estimate for the vanishing LE of a saddle is more difficult than for an attractor.This is because the trajectories considered for the evaluation of the LEs have a finite lifetime, as they unavoidably end up veering towards one of the competing attractors, realizing a positive stretching along the direction of the flow.We have verified (Fig. S1 in the Supplementary Information [SI]) that the estimate of λ 2 converges towards zero as we consider longer-lived trajectories for our estimates, so that the stretching along the flow is minimized (being zero if we could generate a trajectory living exactly on the saddle).
The MLE of our model can be attributed to the chaotic motion induced by the atmosphere via a simple scaling test.First, we note that the magnitude of the maximum and minimum LEs (|λ 1 | and |λ 5 |) is large compared to the other LEs.Furthermore, the values of λ 1 and λ 5 , respectively, are similar for all three invariant sets.This is consistent with the assumption that the atmospheric and oceanic components of our model are weakly coupled, and we may associate λ 1 and λ 5 with the maximal and minimal LE of the uncoupled L84 model, rescaled with the atmospheric timescale 1/ε f .This can be confirmed by computing the Lyapunov spectrum of the attractors for different values of ε f (Fig. S2 in the SI).We find that λ 1 and λ 5 scale in very good agreement with λ i,L84 /ε f , where λ i,L84 is the corresponding LE of the uncoupled L84 model from Eq. 2, while the second to fourth LEs are mostly independent of the timescale separation.
The results shown here exemplarily for σ 0 = 0.9 are qualitatively representative for all values of σ 0 which we have tested within the bistable regime, and there is no clear dependence of λ 1 or λ 2 on the freshwater flux.The MLE λ 1 is very similar across different σ 0 (mean: 536 ± 4) and the second LE is very small (λ 2 < 0.3 ≪ λ 1 ) for all σ 0 .

Dimension of the fractal basin boundary
Given the spectrum of Lyapunov exponents and the lifetime of the saddle, we can use the dimension formula conjectured by [42,43] to relate these two invariants to the information dimension of the stable manifold of the saddle, i.e., the basin boundary.However, provided that the escape rate from the saddle κ = 1/τ is much smaller than the MLE λ 1 , Bódai and Lucarini [44] showed that the formula for the dimension of the basin boundary D b can be simplified to (see their Eq.25) This relation can be viewed as a generalization of the formula proposed earlier by Hsu et al. [73] to the case of rough basin boundaries [44].Note that, strictly speaking, Eq. 6 gives the information dimension and not the box-counting dimension of the basin boundary which will be needed later for assessing the question of predictability, but the two are expected to yield very similar values here [42].
Since the Lyapunov spectra and lifetimes obtained above clearly satisfy κ < λ 1 , we can apply Eq. 6 to obtain a theoretical value for the basin boundary dimension.The resulting dimension D b is between 4.997 and 4.999 for different values of σ 0 , with slightly increasing values as the critical values σ 0,c2 is approached (Fig. S3 in the SI).This means that the basin boundary dimension is extremely close to (within 0.1%) but strictly smaller than the phase space dimension D = 5.Following Eq. 6, this can be attributed to the escape rate being small compared to the MLE.
To verify these results numerically, we also compute the dimension of the basin boundary directly using the standard box-counting algorithm.We follow the methodology of [46] and [44] in sampling evenly spaced initial conditions along a line that intersects with the basin boundary, and determine the attractor to which each initial condition eventually converges.
In Fig. 10a, we show a sample plot of the outcomes of 2 13 initial conditions, where the endpoints of the sampling interval are within δ 1 = 10 −3 of the basin boundary (using the definition of the norm from Sec. 4.1), but away from the saddle.Repeating this procedure in different regions of phase space and then for each value of σ 0 considered previously, we find that the standard box-counting dimension (0 < D box ≤ 1) as shown exemplarily in only increases by 1-2% for a doubling in the accuracy of the initial condition.Note that, formally, this vanishing predictability refers to the full five-dimensional system, but the two competing attractors separate mostly along the oceanic variables (T, S) (cf.Fig. 7).Hence, we can interpret our results as vanishing predictability of the AMOC state.In summary, the final AMOC state is essentially unpredictable in an extended region of phase space close to the basin boundary, to first order irrespective of the value of σ 0 .

Discussion
Using a conceptual climate model that comprises a chaotic atmosphere and a bistable AMOC, we have investigated two fundamental limits to predictability of the asymptotic AMOC state: a fractal basin boundary in the bistable AMOC regime and chaotic transients in the monostable regime.In the bistable regime, we approximated the chaotic saddle (or Melancholia state) between the two competing AMOC attractors using the edge tracking algorithm and found a large timescale separation between the fast chaotic motion on the saddle and the slow escape rate from the saddle.This timescale separation implies a fractal basin boundary with close to full phase space dimension, which we verified by explicitly computing the box-counting dimension of the boundary.In the monostable regime, chaotic transients with exponentially distributed lifetimes arise close to the bifurcation points over a relatively wide range of the freshwater parameter σ 0 .
This complex behavior complicates the assessment of a system's resilience in two distinct ways.In dynamical systems theory [74], one common definition of resilience relies on measuring the minimal "kick" perturbation that causes the system to transition into a competing (undesired) attractor [75].A fractal basin boundary, as quantified by Eq. 1, thus implies a complete loss of resilience in an extended region in state space, where an arbitrarily small perturbation may cause a critical transition.In parameter space, resilience is often defined via the distance to critical points delimiting a "safe operating space" [21,74].The presence of long chaotic transients on "ghost attractors" renders it essentially impossible to determine the exact position of such bifurcation points on finite timescales by observing or simulating only a single time series.Hence, defining a safe operating space requires a probabilistic def-inition due to the exponential distribution of transient lifetimes, and predictability depends crucially on the timescale of interest.
The two phenomena described here have recently been framed as fractality-induced tipping and transient-reduced tipping, respectively, by Kaszás et al. [76].To our knowledge, this is the first study in which both phenomena have been systematically explored for a deterministic but chaotic (conceptual) climate model.In contrast to [76], we did not focus on individual tipping probabilities, but rather on a global characterization of phase space through the uncertainty exponent α (in the bistable regime) and the transient lifetime ⟨τ ⟩ (in the monostable regime).In the bistable regime, this can be achieved with a much smaller computational cost than the local phase space sampling approach taken by [76].In fact, calculating α via the saddle properties and Eq. 6 also appears to yield more accurate results than repeatedly sampling the box-counting dimension of the basin boundary.While both global and local approaches would ideally complement each other, the construction of the saddle and computation of its properties might be computationally feasible in GCMs [46], while a sufficiently large ensemble might not.
To our knowledge, this is the first time that the edge tracking algorithm has been applied to find the chaotic saddle between two competing AMOC states.While it is not guaranteed that our results from a conceptual climate model will generalize to (much) higher-dimensional atmosphere-ocean models, we note several similarities with the results of Lucarini & Bódai [46], who used an intermediate-complexity climate model with O(10 4 ) degrees of freedom to study a different tipping point.From a technical standpoint, the edge tracking algorithm worked robustly, independently of the exact nature of the positive feedback and regardless of whether internal variability of the target variable was larger (as in our case) or smaller than δ 1 .This makes us confident that the edge tracking algorithm is a suitable method even for models with many more degrees of freedom in which bistability of the AMOC has been identified (e.g., [5,25,77,78]), even to those featuring large internal variability [78].Indeed, we have recently implemented the edge tracking algorithm successfully in a coupled atmosphere-ocean GCM of intermediate complexity to investigate a Melancholia state separating the two competing AMOC states [79].While other methods such as numerical continuation have previously been applied to find unstable AMOC equilibria [80], the edge tracking method only requires (sufficiently long) forward integration, which makes it easy to apply and very suitable for tipping elements in explicit models like typical GCMs.
A perhaps surprising similarity to [46] is that the basin boundary dimension is similarly close to full dimension despite the very different models and feedbacks.We believe that this is indeed the case for a rather wide class of models that feature (explicit or emergent) timescale separation, with the local instability governed by the fast component (atmosphere/weather) and the global instability governed by the slow component (ocean/climate).Thus, it seems plausible that sensitive and seemingly random dependence on the initial condition of the final AMOC state [25,81] or of the transient lifetime of a weak AMOC state [82] may be linked to the presence of a high-dimensional chaotic saddle and a fractal basin boundary.
So far, we have only investigated the autonomous case, which is approached when a parameter (here, the freshwater flux σ 0 ) is varied very slowly compared to the internal timescales of the system.This assumption is, however, invalid for the current anthropogenic warming and associated changes in the water cycle.Therefore, rate-dependent AMOC tipping could be observed, but the inertia of the system can also give rise to "safe overshoots" [66] beyond the effect of transient phenomena discussed here.It remains an intriguing open question how the barriers to final state predictability discussed here would be altered in a non-autonomous setting and how they would depend on the forcing rate.

Conclusions
Our results show that predictability of the asymptotic AMOC state is limited by a fractal basin boundary with almost full phase space dimension as well as long chaotic transients, which both arise from the chaotic multiscale nature of the coupled atmosphere-ocean system.
We derived these findings from a conceptual climate model that consisted of four main ingredients for this behavior -bistability, chaotic motion, timescale separation and weak coupling -, such that our conclusions should generalize to more complex models with similar properties and even other components of the Earth system.This has practical implications for the boundaries of a safe operating space, which become intrinsically fuzzy both in a spatial and temporal sense.Starting from a state in the "gray zone" of phase space near a fractal basin boundary or near the "ghost attractor" in the parameter region of chaotic transients, the final outcome depends sensitively on the initial condition, which may give rise to non-monotonic and seemingly counter-intuitive outcomes of an initial condition ensemble -interpretable, however, via dynamical systems theory.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.As discussed in Section 2, ε f is a small parameter that controls the timescale separation between the dynamics of the atmospheric and oceanic component of the system.It is possible to cast Eqs. 4 in the following form, with ϵ 2 = ε f , where x = (x, y, z) and X = (T, S) are vectors comprising the fast and slow variables, respectively.Because of the non-vanishing term g 2 , the system is not skewsymmetric.Hence, one cannot straightforwardly apply the homogenization theory presented in [60,61] to derive the effective stochastic differential equation (SDE) in Itô convention, dX = ã(X)dt + σ(X)dW , (B.2) describing the properties of the slow variables X in the ϵ → 0 limit.Here ã(X) ∈ R 2 is the drift term, σ(X) ∈ R 2×p is the matrix describing the noise law, and dW ∈ R 2 is a vector whose components are the increments of 2 independent Wiener processes.
Instead, one needs to resort to a more general theory [59] that is able to deal with the conceptually more challenging case of two-way coupled systems for rather general functions g(x, X), a(x, X), and b(x, X) in Eq.B.1 above.We derive the following explicit expressions for the drift term and the noise law in Eq.B. where ⊗ indicates the tensor product, µ X (x) is the invariant measure of the system ẋ = g 1 (x) + g 2 (X), where X is taken as a fixed parameter, and x(t) is the evolution at time t of the initial condition x(t = 0) = x.Hence, the correlation matrix of the noise is built starting from the properties of correlation of the fast system where the slow variables X are frozen in.While the general homogenization formulas presented in [59] are more complex, the rather simple result given in Eq.B.3 comes from the fact that in our case the dependence on the fast and slow variables factors out in the functions g, a, and b.Indeed, one obtains the same result one would derive by naively using the theory developed for skew-symmetric systems, under the heuristic assumption that the slow variables can be treated as external parameters for the stochastic process associated with the fast variables.
We also wish to point out the fundamental difference between the original system given in Eq. 4 and its homogenized version for the slow variables given here in Eq.B.2.In the parametric range where multistability is found for (4), it is by definition impossible for a trajectory initialized in the basin of attraction of one state to visit the competing state.By contrast, the stochastic differential equation given in Eq.B.2 allows for (typically rare, when far from the bifurcation points) noise-induced transitions between the competing states of the slow variables (see discussion in [83,84,85]).
Supplementary material for: Limits to predictability of the asymptotic state of the Atlantic Meridional Overturning Circulation in a conceptual climate model Figure S3: Basin boundary dimension calculated from the saddle lifetime and its MLE following Eq.6.
Error bars take into account the uncertainty estimates of the MLEs and the lifetimes.Note that all values are within 1‰ of the full phase space dimension D = 5.
Figure S4: Same as Fig. 10a, but for the full ensemble of 2 13 initial conditions.Zoom in to view the individual lines.
Figure S5: Fraction of initial conditions coverging to the "on" state (black line, blue shading) and to the "off" state (grey shading) along the phase space section shown in Fig. S4.A Lanczos filter with a bandwidth of 150 neighboring initial conditions is used for smoothing.The fraction of initial conditions converging towards the "On" state is clearly, though not monotonously, decreasing as one moves along the phase space section.

Figure 1 :
Figure 1: Bifurcation diagram of the uncoupled Stommel model (Eq.3): fixed points of the AMOC strength Ψ as a function of the control parameter σ.The saddle-node bifurcation points are labeled as L 1 and L 2 .

Figure 2 :
Figure 2: Sketch of the coupled L84-Stommel model.In the Stommel model, T e and S e are temperature and salinity in the subtropical North Atlantic box, and T p and S p are temperature and salinity in the subpolar North Atlantic box, with T ≡ T e − T p and S ≡ S e − S p .The atmosphere-ocean coupling terms via heat fluxes and freshwater fluxes are sketched in red and blue, respectively.

Figure 3 :
Figure 3: An example of long transients for σ 0 = 0.932.Initial conditions were sampled from the "AMOC on" attractor for σ 0 = 0.926 and were integrated forward after the parameter change.The lifetime of the transients spans one order of magnitude and appears unpredictable, but eventually all trajectories convergence to the only remaining attractor.

Figure 4 :
Figure 4: Numerically estimated mean lifetime ⟨τ ⟩ of chaotic transients, corresponding 95% confidence intervals and power-law fit for different values of σ 0 for (a) the transient "AMOC off" state and (b) the transient "AMOC on" state outside the regime of bistability.Note that all axes are logarithmic and that the subplots use different x-axes.

Fig. 4
Fig. 4 shows the lifetime of the transients tracking the "ghost attractor" as a function of the distance to the critical value |σ 0 − σ 0,c |.Here, we only evaluate the range of σ 0 in which the lifetimes of an ensemble of n = 100 initial conditions approximately follow an exponential distribution, and in which τ < T max = 10 4 for all trajectories.Using the fitted critical value, there is a clear power-law relation (a straight line in the log-log plot) between

Figure 5 :
Figure 5: Probability of survival in the region of chaotic transients beyond the critical value σ 0,c2 ≈ 0.927 as a function of the distance to the critical value |σ 0 − σ 0,c | and the timescale of interest (reference time) during which tipping should be prevented.The red marker indicates σ 0 ≈ 0.934 for which the transient tipping probability is limited to 5% on a timescale of 100 years.

Figure 6 :
Figure 6: Schematic depiction of the edge tracking algorithm.Note that in this study, the distances δ 1,2 are calculated over the entire phase space instead of only Ψ. Adapted from Lucarini and Bódai [46], released under a Creative Commons license (CC-BY 3.0, https://creativecommons.org/licenses/by/3.0/).

Figure 7 :
Figure 7: Three-dimensional projection of a pseudo-trajectory of the saddle (green) and trajectories tracking the attractors (blue, orange) for σ 0 = 0.9.Axis show the slow variables T and S and the atmospheric energy

Figure 8 :
Figure 8: Bifurcation diagram for AMOC strength Ψ of the coupled model.Each dot represents one (pseudo-)trajectory averaged over about 50 time units after a spin-up period.Points in the transient regime (gray background) as determined in Sec. 3 are marked with unfilled circles.Vertical solid and dashed lines represent best estimates ± one standard deviation of the critical values σ 0,ci .The colored, solid lines represent square-root fits for the upper and the unstable branches and a linear fit for the lower branch.The bifurcation diagram of the uncoupled Stommel model (Eq.3) is shown in gray for comparison.Saddle-node bifurcation points for the Stommel model are labeled L 1 and L 2 .

Figure 9 :
Figure 9: Lifetimes of the chaotic saddle.(a) Fraction of trajectories N (t)/N 0 within the saddle bounding box B as a function of time for different values of σ 0 .Colored lines indicate linear fits in lin-log space, whose slope is −κ, in the range N/N 0 ∈ [0.03, 0.7].(b) Lifetime τ = 1/κ obtained from the fits in panel a for all simulated values of σ 0 .Error bars reflect the uncertainty in τ with respect to reasonable changes in the fitting range in panel a.The vertical lines and shaded areas indicate estimates of the critical values σ 0,c plus and minus one standard deviation, corresponding to the vertical solid and dashed lines in Fig. 8.

Fig. 9a
Fig.9ashows N S (t)/N 0 for several values of the freshwater forcing σ 0 in the bistable regime of the Stommel model.As expected, the number of remaining trajectories in B shows a clear exponential decay.The lifetime is of O(1) (about one century in physical units) but depends strongly on σ 0 : approaching the saddle-node bifurcation of the underlying Stommel model, trajectories spend more time in the vicinity of the saddle.While the number of sampled σ 0 is too small to conclusively establish a functional relation ⟨τ ⟩(σ 0 ), the lifetime increases clearly non-linearly as the critical freshwater forcing parameter σ 0,c2 is approached (Fig.9b).

Fig. 10b is− 1 2 αFigure 10 :
Fig. 10b is very close to one (between 0.988 and 0.999) for all σ 0 .Both methods therefore confirm, regardless of whether the information or box-counting dimension is used, that a large scale separation between κ and λ 1 yields a basin boundary with almost full phase space dimension.The basin boundary dimension D b and the uncertainty exponent α that determines the final state sensitivity as given by Eq. 1 are linked by the simple relation α = D − D b [40].With the values of D b predicted from the saddle properties (Fig. S3 in the SI), α varies between 0.003 and about 0.001 close to the critical value σ 0,c2 .From the viewpoint of Eq. 1, this means that an improvement in accuracy of the initial condition has a negligible effect on the uncertainty of the final climate state: if the accuracy of the initial condition is doubled, the phase space volume in which the outcome is uncertain decreases by about 1 − 1 2 α ≈ 0.1-0.2%.Even if slightly larger values of α obtained from the box-counting algorithm are used, which are often on the order of 10 −2 , predictability of the final state