# Melancholia States of the Atlantic Meridional Overturning Circulation

Johannes Lohmann<EMAIL_ADDRESS>Physics of Ice, Climate and
Earth, Niels Bohr Institute, University of Copenhagen, Denmark  Valerio
Lucarini  School of Computing and Mathematical Sciences, University of
Leicester, Leicester, UK

###### Abstract

The Atlantic Meridional Overturning Circulation (AMOC) is a much studied
component of the climate system, because its suspected multistability leads to
tipping behaviour with large regional and global climatic impacts. In this
paper we investigate the global stability properties of the system using an
ocean general circulation model. We construct an unstable AMOC state, i.e., an
unstable solution of the flow that resides between the stable regimes of a
vigorous and collapsed AMOC. Such a solution, also known as a Melancholia or
edge state, is a dynamical saddle embedded in the boundary separating the
competing basins of attraction. It is physically relevant since it lies on the
most probable path of a noise-induced transition between the two stable
regimes, and because tipping occurs when one of the attractors and the
Melancholia state collide. Its properties may thus give hints towards physical
mechanisms and predictability of the critical transition. We find that while
the AMOC Melancholia state as viewed from its upper ocean properties lies
between the vigorous and collapsed regimes, it is characterized by an Atlantic
deep ocean that is fresher and colder compared to both stable regimes. The
Melancholia state has higher dynamic enthalpy than either stable state,
representing a state of higher potential energy that is in good agreement with
the dynamical landscape view on metastability.

##  I Introduction

The Earth’s climate history features periods of smooth changes interlaced with
rapid variations likely due to critical transitions (Westerhold  _et al._ ,
2020; Ghil and Lucarini, 2020; Boers  _et al._ , 2022; Rousseau  _et al._ ,
2023), the latter being sometimes associated with mass extinctions events and
fundamental changes in the biosphere (Lenton, 2013; Alsulami and Petrovskii,
2023; Algeo and Shen, 2024). Critical transitions in the Earth system arise
when the states of different large-scale climatic sub-components lose
stability. There is widespread concern that ongoing climate change could
indeed trigger the onset of transitions in various so-called tipping elements
of the climate system. These include the Amazon Forest, the Atlantic
Meridional Overturning Circulation (AMOC), the Greenland ice sheet, and the
Permafrost, among others (Lenton  _et al._ , 2008; Armstrong McKay  _et al._ ,
2022). Hence it is nowadays rather common to refer to the ongoing impact of
anthropogenic activities on the climate system as the climate crisis (Lucarini
and Chekroun, 2023).

The Earth system consists of a variety of sub-systems that evolve on time
scales spanning many orders of magnitude (Ghil and Lucarini, 2020; von der
Heydt  _et al._ , 2011). The sub-systems are coupled by various physical
processes, and as a result the combined Earth system is a multi-scale
dynamical system. But it has long been argued, along the lines of the
Hasselmann program (Imkeller and von Storch, 2001), that in some circumstances
these time scales may be treated as sufficiently separated, so that a sub-
system evolving on a particular time scale can be considered to be forced by
stochastic variability in sub-systems that evolve on much faster time scales
(Hasselmann, 1976). This motivates why stochastic parametrizations are needed
to improve the performance of weather forecast and climate models (Berner  _et
al._ , 2017). One the other extreme, if there are sub-systems that evolve on
even much longer time scales than the element in question, the effect of these
sub-systems may be considered a slowly-varying, quasi-adiabatically changing
control parameter on the time scale of interest (Saltzman, 2001). If one
relaxes the assumption of time-scale separation, the Hasselmann program
requires a careful re-examination, because rates of change become non-
negligible, the stochastic forcing is not white, and the dynamics includes
non-markovian contributions (Lucarini and Chekroun, 2023).

As clarified in Ashwin  _et al._ (2012), tipping behaviour can result, by and
large, from three different mathematical mechanisms, which are often co-
existent in real-life systems: 1) classical bifurcations associated with
quasi-adiabatic modulations of the parameters of the system (B-tipping); 2)
transitions occurring specifically as a result of non-negligible rates of
change of one or more system parameters (R-tipping); and 3) stochastic
perturbations inducing transitions between competing modes of operation of the
system (N-tipping).

These tipping mechanisms all rely on an underlying multistability of the
system, which originates in nonlinear dynamical laws associated with positive
feedbacks, leading to the existence of multiple competing stable states. A
notable example of multistability in the Earth system is the AMOC, the
stability properties of which have long been studied beginning with the
landmark contribution by Stommel (Stommel, 1961), who highlighted using a
simple box model that the competing (destabilizing) salinity and (stabilizing)
temperature advection feedbacks can result in more than one competing state.
Here, the freshwater flux into the ocean is the natural parameter controlling
whether the system is in a monostable or bistable configuration. One state -
resembling the present-day AMOC conditions - is characterized by a vigorous
interhemispheric circulation (’ON’ state). The competing state instead
features weak or absent large scale flow (’OFF’ state). The driving
mechanisms, controls, and stability of the AMOC are being actively studied
(Kuhlbrodt  _et al._ , 2007; Weijer  _et al._ , 2019), especially because
variations in the intensity and features of the AMOC can have key impact on
the regional and even global climate (Rahmstorf, 2002; Liu  _et al._ , 2020),
and are thought to have played a major role in large climatic variations in
the Pleistocene (Broecker  _et al._ , 1990; Barker  _et al._ , 2021). A wealth
of studies performed using models of different complexity have confirmed the
overall picture proposed by Stommel (Rahmstorf  _et al._ , 2005; Hawkins  _et
al._ , 2011), and observations of so-called early-warning signals associated
with B-tipping indicate that the ongoing climate change might be leading to a
loss of stability of the current active state of the AMOC, i.e., a loss of the
present intense interhemispheric ocean heat transport (Boers, 2021; Ditlevsen
and Ditlevsen, 2023; van Westen  _et al._ , 2024).

A recent modelling study (Lohmann  _et al._ , 2024) indicates that the AMOC
might possess, for certain values of the system parameters, well more than
just two competing states: several distinct variants are found within both
regimes of strong and weak circulation, associated with a rather complex
bifurcation diagram and a web of possible transitions between the various
competing states as parameters are varied. The presence of such a complex
dynamical landscape suggests that the above-mentioned early-warning signals do
not necessarily indicate that the loss of stability of the current state will
lead to the AMOC falling into the ’OFF’ state. While early-warning signals
investigate local properties of the system, it is clear that a more global
view on the stability landscape of the system is helpful.

###  I.1 Multistability and Edge States

In multistable systems, the asymptotic state depends on the initial
conditions, with different subsets of initial conditions (basins of
attraction) evolving towards different stable states (attractors). If the
system is initialised on the boundary between neighbouring basins of
attraction, the state evolves towards a special asymptotic set, the edge
state, which is a saddle that separates the competing attractors. Small
perturbations away from the saddle result in the system converging to either
attractor with probability one.

Edge states in physical systems were first studied in the context of classical
fluid dynamical configurations like the plane Couette flow and the pipe flow
(Schneider  _et al._ , 2008; Schneider and Eckhardt, 2009), even if,
rigorously speaking, such systems possess just one true attractor (the laminar
flow), which is co-existing, in the high-forcing regime, with an extremely
long-lived turbulent transient (Hof  _et al._ , 2006, 2008). The so-called
edge tracking algorithm (Skufca  _et al._ , 2006; Vollmer  _et al._ , 2009)
has been a key tool to construct edge states, which are otherwise inaccessible
using standard forward numerical integration because of its unstable nature.

The edge tracking algorithm has also allowed the computation of edge states in
nontrivial climate models featuring highly turbulent behaviour and co-existing
warm and snowball states (Lucarini and Bódai, 2017). There, the dynamical
saddle was found to be a chaotic state, and in this particular case was
referred to as Melancholia (M) state, a nomenclature we will follow hereafter
as the dynamics in the system we consider here is also chaotic. The presence
of a chaotic saddle was also shown to be linked to the presence of a fractal
basin boundary between the competing states, featuring a co-dimension strictly
smaller than one (Vollmer  _et al._ , 2009; Bódai and Lucarini, 2020). As a
result, one finds a vast region of phase space where the so-called
predictability of the second kind (Lorenz, 1975) is extremely limited:
applying arbitrarily small perturbations to the initial conditions can lead to
system ending up in a different asymptotic state (Mehling  _et al._ , 2024).

The M state is physically relevant because tipping behaviour associated with
quasi-adiabatic parametric modulations occur when one of the attractors and
the M state collide. Hence, its properties may give hints towards physical
mechanisms and predictability of the critical transition. The M state may also
be key to interpret some cases of rate-induced tipping, where a temporary
attraction to the M state can occur (Wieczorek  _et al._ , 2011), resulting in
long transients (Lohmann  _et al._ , 2021; Lohmann and Ditlevsen, 2021).

Finally, the M state can be of great relevance in the case of noise-induced
tipping. The presence of stochastic forcing allows for transitions between the
different basin of attraction, leading to the phenomenon of metastability.
Graham’s field theory (Graham, 1987; Graham  _et al._ , 1991) indicates that
in the weak-noise limit transitions between competing basins of attraction can
be described - both in terms of time statistics and paths - by large deviation
theory (Touchette, 2009). Transition paths are concentrated around a special
trajectory known as instanton, which connects the attractors of the system
with the M state living on the basin boundary. The M state can thus be
considered the gateway of noise-induced transitions. Identifying the physical
characteristics of an M state or instanton can be useful as an indicator that
a transition is on the way, and it can help identify the physical mechanisms
and perturbations that promote a transition.

###  I.2 This paper

Evidence for regime transitions in geophysical flows that may be understood as
stochastic transitions has been found for oceanic western boundary currents
(Schmeits and Dijkstra, 2001), zonal-blocking transitions in the atmosphere
(Plaut and Vautard, 1993), and geomagnetic reversals (Glatzmaier and Roberts,
1995; Valet, 2003). Previous work has identified instantons and edge states
from simulations of geophysical flows, such as for atmospheric jet dynamics
(Bouchet  _et al._ , 2019) and the snowball Earth transition (Lucarini and
Bódai, 2019, 2020; Margazoglou  _et al._ , 2021). But the M state associated
with regime transitions of the AMOC has only been started to be investigated
recently (Börner  _et al._ , 2024). A case for time scale separation,
stochastic forcing and spontaneous transitions of the AMOC has been made
previously (Cessi, 1994; Drijfhout  _et al._ , 2013; Kleppin  _et al._ , 2015;
Castellana and A., 2020), and the possibility of a spontaneous AMOC collapse
has been recently observed using a rare-event algorithm in the case where the
turbulent flow of the atmosphere was explicitly simulated in a coupled climate
model (Cini  _et al._ , 2024). In such a case, the AMOC is not, strictly
speaking, multistable, because the system performs a transition from the state
with vigorous AMOC to the state featuring weak AMOC in the absence of any
stochastic forcing and purely as a result of internal variability.
Nonetheless, the residence time of the system in the vigorous AMOC state is
extremely long.

Here we compute the M state that lies between the vigorous and collapsed AMOC
regimes, based on the recently constructed stability landscape of a global
ocean model under North Atlantic (NA) freshwater forcing (Lohmann  _et al._ ,
2024). Its physical characteristics are analyzed in order to assess whether
they can be harnessed to detect if a spontaneous collapse of the AMOC is
underway, and to elucidate the physical mechanisms that would facilitate such
a transition.

The structure of the paper is as follows. First, we present a qualitative
survey of the co-existing attractors and their basins of attraction based on a
large ensemble simulation (Sec. II.1). Next, simulations are performed that
precisely locate the basin boundary separating the vigorous and collapsed AMOC
regimes, and, from a given position on this basin boundary, edge tracking is
performed (Sec. II.2). The physical properties of the computed M state are
analyzed in Sec. II.3, and conclusions are drawn in Sec. III. Finally,
Appendix A contains details on the numerical model used here and on the
numerical simulations performed in this study.

##  II Results

Figure 1:  (a) Stability diagram of AMOC in the ocean model, as obtained in
Lohmann  _et al._ (2024). Symbols of different color indicate the different
branches of attractors. Dotted lines and their color indicate which attractor
branch collapses onto which other branch. The vertical line indicates the
parameter value at which we performed ensemble simulations that are presented
in this study. (b) Maximum AMOC strength time series of the individual
realizations of the ensemble (gray) with $N=199$ members, along with longer
steady-state simulations (at least 10 000 years) started close to each of the
co-exisiting attractors (colored time series).

###  II.1 Survey of the basin structure

The AMOC as simulated by the primitive equation finite-difference global ocean
model Veros (Häfner  _et al._ , 2018, 2021) has been shown to display a high
degree of multistability (Lohmann  _et al._ , 2024), as seen in the
bifurcation diagram in Fig. 1a, where symbols of same shape and color indicate
the different branches of attractors and the control parameter on the
horizontal axis is the NA freshwater forcing $F$ (see A). We focus here on the
forcing parameter $F=F_{0}=0.347$ Sv, which lies in the middle of the main
hysteresis loop (vertical line in Fig. 1a). At $F_{0}$ there is a cluster of
three distinct states featuring a very weak AMOC, as well as a cluster of
states featuring a vigorous AMOC. This second cluster is further divided into
two bands of states with an average strength of approximately 8.2 and 8.7 Sv,
respectively (see colored time series in Fig. 1b and inset). For the band with
slightly larger AMOC strength, three distinct attractors have been identified
in Lohmann  _et al._ (2024) (see inset Fig. 1a). In the lower band only one
attractor was identified, but both the upper and lower bands might contain
more attractors. For slightly higher $F$, an additional cluster of three
partially-collapsed states forms another smaller-scale hysteresis loop (inset
in Fig. 1a), in agreement with the hierarchical structure of the dynamical
landscape of the climate system conjectured in (Margazoglou  _et al._ , 2021).
In terms of different macroscopic observables, states within the clusters may
be separated more clearly. The open circles in Fig. 2a show a two-dimensional
projection of the mean states on the attractors, where the cluster of
collapsed AMOC states is seen to split up in terms of the NA subsurface
temperature. Similarly, the two bands of vigorous states are well-separated in
terms of other observables (see e.g. open symbols in Fig. 2c).

As shown in Lucarini and Bódai (2017) and Mehling  _et al._ (2024), computing
the M state of a bistable system is relatively straightforward once one has a
working edge tracking algorithm in place, the main challenge being the need to
initialize the procedure with two very nearby initial conditions on the two
opposite sides of the basin boundary. Against intuition, the procedure works
even in the case of M states, where the basin boundary has codimension
strictly smaller than one - and possibly very close to zero. In the case
studied here, the situation is considerably more complicated because the
system has at least seven distinct chaotic attractors, which leaves the door
open for a large variety of geometrical configurations for the basins of
attraction and of the features of the basin boundaries.

Hence, as preliminary step to computing an M state, we perform a survey of the
basins of attraction and their boundaries for the co-existing attractors,
using ensemble simulations with a large variety of initial conditions. This is
achieved by chosing states along several straight lines in phase space,
yielding a handful of one-dimensional ’sections’ of phase space. The lines we
construct are defined by two states $A$ and $B$ on a pair of attractors, where
one attractor features a vigorous AMOC and the other a collapsed AMOC. The
states are chosen as the model snapshot at the last time step of the long
steady-state simulations that were performed at fixed forcing values in
Lohmann  _et al._ (2024), and thus approximately correspond to one individual
point in each of the distinct attractors. We then perform a linear
interpolation of the entire set of state variables and construct 32 equally
spaced initial conditions $C$:

$C=\alpha A+(1-\alpha)B,\,\,\,\,\text{where }\alpha=(i-5)/20,\,\,\,\,\text{and
}i=\\{0,1,..,31\\}.$ (1)

The values $i=\\{5,6,..,25\\}$ correspond to an interpolation between the pair
of states, and considering larger and smaller values of $i$ allows us to
sample a larger variety of initial ocean states.

Figure 2:  (a) Two-dimensional projection of the asymptotic states of the
ensemble members onto the maximum AMOC strength and the NA sub-surface
temperature, averaged over the last 400 years of the simulations (black dots).
The open, colored circles are 2000-year mean values on the co-existing
attractors. We labeled the five of the attractors (I to V) that have been used
to construct the ensemble simulations. (b) Sequence of basins of attraction
along six different lines L1-L6 across phase space, defined by interpolating
between pairs of states on the attractors II-V, I-V, I-III, II-III, II-IV, and
I-IV, respectively. (c,d) Initial states of the ensemble members projected
onto two scalar observables, and color coded according to which attractor the
simulations converge to (for color coding see Fig. 1). Projection onto maximum
AMOC strength and South Atlantic subsurface salinity (c) and temperature (d).

To construct initial conditions according to Eq. 1, we choose the three weak
AMOC attractors and two of the four vigorous AMOC attractors (dark blue and
red time series in Fig. 1b). We are unable to consider all four vigorous
attractors (cyan and light yellow are excluded) because of computing
limitations. Nonetheless, we are able to cover both bands of vigorous
attractors, see above. This yields six pairings of a vigorous and a collapsed
attractor and an ensemble comprising initial conditions lying on six different
lines L1-L6 in phase space (see Fig. 2a,b for definition of the attractors and
lines). Figure 1b shows how the ensemble trajectories (gray) converge to the
different attractors (colored time series). Within the simulation time of 3
000 years it can be determined from the AMOC strength which attractors each
ensemble member converges to, with the exception of not being able to clearly
distinguish the close-by states in the uppermost band of the vigorous AMOC
attractors (inset Fig. 1b).

Figures 2a show a projection of the asymptotic states of the ensemble members
(black circles) defined by the mean AMOC strength and NA sub-surface
temperature. The ensemble converges to the two distinct bands of vigorous AMOC
attractors, but within the bands a distinction of attractors is difficult
(inset Fig. 2a). Longer simulation times would be required to converge even
closer to the attractors and to allow for longer time averages that yield more
precise mean values by eliminating the centennial-scale variability. As a
result, we neglect some of the fine structure and in the following treat the
uppermost band as one attractor color coded in dark blue.

Pairing the initial conditions of the ensemble members with their asymptotic
states gives useful information on the basin boundaries for $F=F_{0}$. Figure
2c,d shows two-dimensional projections of the initial conditions that are
color coded according to which attractor they converge to. The sequence of
basins along the individual lines in phase space that constitute the ensemble
is shown as a function of $\alpha$ (Eq. 1) in Fig. 2b. The connection from
vigorous to collapsed AMOC regimes is always via the sequence dark blue basin
- green basin. These are the basins of the two attractors we refer to as ’ON’
and ’OFF’ hereafter. Although based on a limited number of initial conditions
and very sparse exploration of an extremely high-dimensional phase space, this
may indicate that as a result of how the basins are geometrically arranged in
phase space, there is a specific path that needs to be traversed in a
transition between the vigorous and collapsed AMOC regimes, regardless of
which specific vigorous and collapsed attractors are the starting and end
points.

We remark that the dark blue and green states are not those one would naively
choose by looking at the AMOC intensity as indicator of proximity: the dark
blue attractor has a more vigorous circulation than the red one; and the AMOC
of the green attractor is slightly weaker than that of the yellow attractor.
This further indicates that the AMOC intensity alone does not necessarily
provide sufficient information to understand the detailed dynamics of such a
high-dimensional flow.

###  II.2 Basin Boundaries and Edge Tracking

To precisely locate the basin boundary, we follow the initial step of the edge
tracking algorithm and begin with an arbitrary pair of initial conditions that
converge to ON and OFF, respectively. We bisect the two initial conditions and
use the resulting state as initial condition for a new simulation, where it is
determined whether the trajectory converges to ON or OFF. Next, we bisect the
newly constructed initial state and the previous initial condition that
converged to the opposite attractor, and use the resulting state as initial
condition for the next simulation. This procedure is repeated for $M$
iterations, which yields two very close-by initial conditions on opposite
sides of the basin boundary.

Figure 3:  Repeated bisections of initial conditions are performed, such that
a pair of very close-by initial conditions is obtained where one reaches the
ON attractor, and the other the OFF attractor. (a) For each iteration of the
bisection procedure, we estimate the time until the pair of trajectories has
separated significantly. This is defined by the first time the trajectories
deviate in AMOC strength by more than 0.2 Sv, as estimated from the time
series smoothened by a 20-year Gaussian kernel filter. (b) Time series of the
AMOC strength of the individual realizations during the bi-section procedure,
drawn with increasing contrast for increasing iteration number. The inset is a
close-up on the time series after the fifth iteration. (c) Two-dimensional
projection onto the AMOC strength and NA sub-surface temperature of the pair
of trajectories at the final iteration.

Figure 3 shows one instance of this procedure starting from a pair of initial
conditions on the line L3 that belong to the dark blue and green basins,
respectively, as defined in Fig. 2. At each iteration, we measure the time
elapsed until the pair of trajectories separates appreciably in terms of AMOC
strength (Fig. 3a), after which they converge to opposite attractors. This
time increases rapidly for the first few iterations and then stabilizes at
several hundred years (Fig. 3b). The final pair of trajectories is very close-
by at $\alpha=\frac{6453063}{16777216}=0.384632528$ and
$\alpha=\frac{806633}{2097152}=0.384632587$, and evolves closely in phase
space for almost 500 years before separating (Fig. 3c). This indicates that
both trajectories are closely following the stable manifold of the M state
(the basin boundary) before the transversal instability takes over (Lucarini
and Bódai, 2017). These trajectories are then used to initialize the main
component of the edge tracking procedure.

We repeatedly apply the bisection method above to new initial states that are
created by a relatively short time evolution of the two closest initial
conditions on opposite sides of the boundary. At each iteration $i$, pairs of
initial conditions that lie close-by on opposite sides of the basin boundary
are constructed by a fixed number of $k$ bisections, in the following way.
Starting from a pair of initial states $\mathcal{A}$ and $\mathcal{B}$ that
converge to the ON and OFF attractor, respectively, a bisection is performed
to obtain a new initial state $\mathcal{C}=(\mathcal{A}+\mathcal{B})/2$.
$\mathcal{C}$ is then evolved from time $t_{0}$ to $t_{2}$, chosen such that
we can decide with certainty from a threshold on the AMOC strength whether the
trajectory will converge to ON or OFF. If $\mathcal{C}$ converges to ON we
assign $\mathcal{C}\rightarrow\mathcal{A}$, and if it converges to OFF we set
$\mathcal{C}\rightarrow\mathcal{B}$. Thereafter a new initial condition
$\mathcal{C}=(\mathcal{A}+\mathcal{B})/2$ is constructed, and the procedure is
repeated until $k$ bisections are reached. The trajectories associated with
the states $\mathcal{A}$ and $\mathcal{B}$ are saved and a new iteration $i+1$
of the edge tracking procedure is spawned by using the model snapshots
obtained after the time evolution of $\mathcal{A}$ and $\mathcal{B}$ until
$t_{1}<t_{2}$ as starting initial states for a new round of $k$ bisections.
$t_{1}$ is chosen such that the trajectories begin to separate, but are not
too far from the boundary yet.

If the bisections yield initial states close enough to the boundary, at each
iteration the model states move a little closer towards the M state along the
basin boundary, as we are shadowing a true trajectory that is moving along the
attractive, stable manifold of the M state. Figure 4 shows results of the edge
tracking procedure initiated from model snapshots taken shortly after the
separation of the trajectories shown in Fig. 3c. We choose $k=7$ bisections,
since at the given forcing $F_{0}$ this is sufficient to generate pairs of
trajectories that evolve along the boundary for a longer period of time (Fig.
3a). Further, $t_{1}=100$ years and $t_{2}=250$ years are used, and all
trajectories are continued until $t_{3}=500$ years to check whether they
indeed continue on the path to the attractor that was determined at $t_{2}$.
We use $t_{2}<t_{3}$ for computational reasons, such that more realizations
can be simulated in parallel. All time series generated during the edge
tracking are shown in Fig. 4a. The algorithm produces pairs of trajectories
that reliably evolve close to each other for roughly a century (Fig. 4b,c),
after which they quickly separate, allowing us to detect reliably at $t_{2}$
whether they converge to ON or OFF. Here we used trial and error to roughly
match $t_{1}$ to the separation time of the closest pairs of trajectories.

Figure 4:  (a) All simulations performed during one instance of the edge
tracking procedure with $N=60$ iterations, shown as time series of the AMOC
strength. The longer red and black trajectories are the simulations from which
the initial state is chosen. Thereafter, the edge tracking algorithm produces
subsequent generations of trajectories, where for each iteration the colored
time series correspond to the pair of simulations with the closest-by initial
conditions that converge to opposite attractors. The gray trajectories are all
other simulations. (b) Closeup of the time series for three subsequent
iterations. The threshold used to decide whether a simulation is converging to
ON or OFF is shown with the black, dashed line. (c) As in (a), but only the
final pairs of each iteration are shown, and only the portion of the time
series until $t_{1}$ (when a new iteration is branched off) is given in color.

It is apparent from Fig. 4 that the AMOC strength does not provide any useful
indication of the motion along the basin boundary, as by looking at its time
series we are unable to elicit whether the true trajectory we are shadowing
has reached the M state. Figure 5a shows a two-dimensional projection of the
trajectories of the pairs of closest initial conditions for all iterations.
The evolution of the AMOC strength is very similar for all iterations, but
other observables related to the deep ocean water properties show a clear
trend as the edge tracking procedure progresses. During early iterations,
subsequent trajectories in Fig. 5a are spaced relatively far apart along the
horizontal axis, as they move quickly towards initial conditions with lower
Atlantic deep ocean salinity. For later iterations subsequent trajectories get
closer and closer together, until they converge to a stationary region in the
two-dimensional projection.

Figure 5:  (a) Pairs of trajectories in 2D projection of AMOC strength and
tropical Atlantic deep ocean salinity for all edge tracking iterations. The
iterations are progressing from right to left, which is further indicated by
increase in contrast of the trajectories. The red (black/gray) trajectories
are the simulations that converge to OFF (ON). The average initial states of
the trajectory pairs converge to the estimated M state, shown by the green
triangle. The estimated basin boundary is given by the blue dashed line. (b)
Average state in terms of AMOC strength of the final pair of trajectories for
each iteration, obtained from the average model state of the first 5
simulation years, averaged over the trajectory pairs. (c) Same, but for the NA
deep ocean temperature.

By averaging over the first 5 simulation years of a pair of trajectories at a
given iteration, we estimate the shadowed state on the basin boundary. As the
edge tracking iteration proceeds, such a state moves along the boundary (blue
dashed line in Fig. 5a) and asymptotically tends towards the M state (green
triangle, estimated from the average state of the last 10 iterations), showing
that the M state is indeed a relative attractor for the basin boundary
(Lucarini and Bódai, 2017). The AMOC strength of the shadowed state does not
change much after a few iterations (Fig. 5b). But there is an approximately
exponential relaxation of the slow variables associated with the deep ocean,
converging to a steady state after 30-40 iterations (Fig. 5c). Given
$t_{1}=100$ years, the resulting relaxation time is in good agreement with the
typical multi-millennial time scale until the model equilibrates onto one of
its attractors (Lohmann  _et al._ , 2024).

###  II.3 Characteristics of the Edge State

We now compare the specific M state that was computed, which we will denote as
’EDGE’ hereafter, to the ON and OFF attractors in terms of climatological
averages. To this end we consider the final 10 iterations of the edge tracking
procedure and average the model state over the first 5 years of the final
pairs of simulations. For the attractors we use the mean state obtained during
the final 500 years of the long steady-state simulations. In terms of the NA
circulation, EDGE lies ’between’ ON and OFF featuring an intermediate strength
of the AMOC and NA deep water formation (Fig. 6a-f). EDGE seems relatively
closer to OFF compared to ON, for instance showing a weakening of both NA
gyres that is almost as pronounced as the weakening in OFF (Fig. 6g-i, note
positive anomalies of the barotropic streamfunction due to its negative sign
in the subpolar gyre).

Further, as opposed to ON, both EDGE and OFF lack any substantial inter-
hemispheric transport (Fig. 6a-c). As a result, EDGE also shows strong cooling
and freshening of the NA surface waters (Fig. 7). Especially the temperature
decrease is almost as pronounced as in OFF. This is not surprising since the
AMOC is indeed severely weakened in EDGE, but it nevertheless confirms that
the AMOC slowdown fingerprints - see (Roberts  _et al._ , 2013; Rahmstorf  _et
al._ , 2015; Caesar  _et al._ , 2018; Jackson and Wood, 2020) \- would also be
observed in the case of an AMOC collapse via EDGE.

Figure 6:  (a-c) The Atlantic overturning streamfunction of the M state (b) is
compared to those of the ON and OFF attractors. The contours show the
circulation strength in Sv. (d-f) Same, but for maps of the average winter
mixed layer depth (MLD) in the NA. MLD is defined by the depth until which the
monthly average ocean temperature is within 0.5 K of the sea surface
temperature, indicating the degree of vertical mixing via convection. The
winter MLD is obtained by averaging over the first 4 months of the year. (g-i)
Same, but for the barotropic streamfunction in the NA.

In the upper 1000 m of the global ocean the spatial pattern of the average
temperature, salinity, and density anomalies of EDGE compared to ON are very
similar to, but slightly weaker than those of OFF compared to ON (Fig. 8a-d).
Most significant is a NA cooling and freshening, and a SA warming and increase
in salinity of EDGE and OFF compared to ON. Thus, in terms of the upper ocean
properties, EDGE lies between ON and OFF. Going from ON to either OFF or EDGE
shows clear a decrease in upper ocean density. Thus, as a consequence of an
AMOC slowdown the NA freshening dominates the NA cooling, which supports that
the destabilizing salt advection feedback is stronger than the stabilizing
temperature advection feedback.

The situation in the deep ocean is markedly different, since here the
properties of EDGE are outside of the range of the ON and OFF regimes for a
large part of the global ocean. In particular, EDGE is colder and fresher than
both ON and OFF in the entire Atlantic expect the Nordic seas, as seen by the
same sign of the anomalies in Fig. 9b,e and c,f. Due to the opposite effect of
temperature and salinity, the pattern for the deep Atlantic density is more
fragmented. In parts of the Atlantic freshening dominates over cooling, so
that EDGE is less dense than both ON and OFF (Fig. 9g-i). However, in other
parts of the global ocean the density of EDGE lies between ON and OFF.

Figure 7:  Anomaly in the sea surface temperature when moving from ON to OFF
(a), and from ON to EDGE (b). (c,d) Same, but for sea surface salinity.

Having found a first indication of a characteristic anomaly in EDGE, we now
want to consider how one can use the knowledge of such an anomaly to identify
an ongoing AMOC collapse. Let us assume we are in the ON state and the AMOC is
slowing down. We further assume that from climate modeling and paleoclimatic
observations we have knowledge on the existence and features of an AMOC OFF
state. We do not know whether the observed AMOC slowdown is only a temporary
excursion, or whether it will lead us dangerously close to tipping. This needs
to be assessed by consulting specific observables or fingerprints. Without
specific knowledge of the exact transition path of an AMOC collapse - since it
has not been directly observed before - it would be a priori assumed that the
fingerprints pointing towards a possible transition from ON to OFF are in the
direction of the ON-OFF anomaly, and that an ongoing transition should occur
in the same direction. But from the theory of N-tipping and B-tipping we know
that the relevant direction is instead determined by the position of EDGE with
respect to ON. Hence, better fingerprints are those that are aligned in the
ON-EDGE direction. While for some observables the ON-EDGE and ON-OFF anomalies
point in the same direction, it is not the case for others. It is most
important to identify the fingerprints where ON-EDGE and ON-OFF are not
aligned, since in these cases the exclusive knowledge of the OFF state would
lead one to believe that there is no dangerous excursion - or even an overall
AMOC strengthening - because the system is moving further away from OFF. Note
that if we are in the OFF state, the relevant direction for a transition is
given by the EDGE-OFF anomaly. Hence, good fingerprints for the ON-to-OFF and
OFF-to-ON transition are in general different.

Figure 8:  Differences between the attractors and the M state in the
temperature (a-c), salinity (d-f), and density (g-i) averaged over the upper
1000 m of the ocean.

With this in mind, for the structure of the deep Atlantic in EDGE we find that
the ON-EDGE and ON-OFF anomalies point in the same direction, since OFF lies
between ON and EDGE. Thus, there is a trend towards fresher and colder deep
Atlantic both when expecting a hypothetical ’direct’ transition from ON to OFF
(Fig. 9a,d), and when observing an excursion from ON towards EDGE (Fig. 9b,e).
The converse holds when considering an AMOC resurgence starting from OFF,
because the OFF$\rightarrow$ON and OFF$\rightarrow$EDGE trends differ in sign.
The AMOC resurgence scenario is not relevant for present-day climate, but its
fingerprints may be useful for understanding abrupt changes in paleoclimate,
such as the Dansgaard-Oeschger events. Outside the Atlantic, there are
extended regions where the deep ocean density in ON lies between the densities
in EDGE and OFF. Thus, in contrast to the Atlantic, in these regions the ON-
OFF and ON-EDGE density anomalies point in different directions. Most notably,
the density of the deep Indian ocean and almost the entire Pacific is lower in
ON compared to OFF, but even lower in EDGE (Fig. 9g,h). Here, a transition
from ON via EDGE towards OFF would first see a decrease in density and then an
increase. The initial density decrease is a non-trivial fingerprint of an AMOC
collapse as it would not be expected from the knowledge of the OFF state only.

Figure 9:  Same as Fig. 8, but averaged for the ocean below 1000 m depth.

As a result of the combined effect of temperature and salinity anomalies on
the density, the deep ocean in EDGE is lighter on average compared to both ON
and OFF (Fig. 9h,i), when disregarding the Southern ocean. The lighter waters
at depth raise the center of gravity, and thus EDGE has the highest
gravitational potential energy of the three states. This can be quantified by
computing the dynamic enthalpy Young (2010), which is the effective potential
energy of the ocean and measures the energy released when moving a fluid
parcel of unit mass adiabatically from its current position to the surface
Nycander (2010). Note that since there is stable stratification most of the
time, this energy is negative on average and energy is required to move
parcels to the surface. Compared to the attractors, in EDGE there is less
energy required to move parcels up and thus it represents a higher potential
energy state. Specifically, the vertically averaged dynamic enthalpy shows
that EDGE is in a higher energy state compared to OFF for the entire ocean
except the Nordic seas (Fig. 10b), and in a higher energy state for the
Atlantic compared to ON (Fig. 10a). In the global average, the dynamic
enthalpy in EDGE is higher by 2.09$e^{20}$J and 4.86$e^{19}$J compared to OFF
and ON, respectively.

Figure 10:  Change in dynamic enthalpy (averaged over the entire water column)
when moving from ON to EDGE (a), and from OFF to EDGE (b).

Based on this, we propose that the global stability properties of the AMOC can
be understood in terms of a quasipotential that can be approximated by the
total dynamic enthalpy in the Atlantic (north of 30oS). A one-dimensional
projection of such a quasipotential is sketched in Fig. 11a, where all values
have been offset such that the potential energy in EDGE is 0. The potential in
OFF is higher than ON, consistent with the observed closeness of EDGE and OFF
over most metrics concerning the Atlantic. This indicates that the system is
closer to the loss of stability of the OFF regime for this value of the
control parameter.

In the weak-noise limit, a typical noise-induced transition between ON and OFF
evolves close to the instanton from the departing attractor to EDGE, and then
follows the relaxation trajectory from the M state to the final target
attractor. Relaxation and instantonic trajectories are different in the case
of nonequilibrium systems (Graham, 1987; Graham  _et al._ , 1991; Lucarini and
Bódai, 2020; Margazoglou  _et al._ , 2021). We obtain full relaxation
trajectories starting very close to EDGE by extending the simulation time for
the pairs of trajectories during the 10 last iterations of the edge tracking
algorithm to 4 000 years, and project them in a two-dimensional observable
space (purple dots in Fig. 11b). Computing instantonic trajectories is beyond
the scopes of our study because it requires either running lengthy
stochastically forced simulations and then performing suitable stochastic
averages (Margazoglou  _et al._ , 2021), or solving a complex variational
problem (Grafke and Vanden-Eijnden, 2019). Nonetheless, based on the location
of the M state and the relaxation trajectories, it is clear that when starting
at either attractor potential energy needs to be gained by bringing relatively
fresh water to the deep ocean. This corresponds to a movement in the negative
vertical direction in Fig. 11b, as indicated by the arrows. After reaching
EDGE, there is a relaxation towards the attractors during which potential
energy is released again (movement in positive vertical direction in Fig.
11b).

Figure 11:  (a) Conceptual quasipotential landscape projected onto the AMOC
strength, derived from the differences in Atlantic dynamic enthalpy of the
attractors and the M state. (b) Two-dimensional projection of the mean states
in the attractors (colored dots) and the M state (green triangle), as well as
the basin boundary (blue dashed line) and the relaxation trajectories
(purple). (c) Time average of the total Atlantic dynamic enthalpy in the
attractors and M state (horizontal lines). The blue time series is a pseudo-
trajectory approaching the M state during the edge tracking algorithm,
obtained by splicing together the 100-year paths (given by $t_{1}$) that
correspond to the movement along the basin boundary for subsequent iterations.
Two trajectories that relax from the vicinity of the M state (starting at
$t=6000$ years) to the ON and OFF attractors are also shown.

The Atlantic dynamic enthalpy captures the physics of EDGE as a state of
higher potential energy compared to both ON and OFF. Yet, it is not an exact
Lyapunov function of the system. Trajectories starting at EDGE show the
expected roughly monotonic decrease in the potential, as seen by two
trajectories (purple and gray time series in Fig. 11c) leading to the
competing attractors. But also we find initial conditions where the
deterministic evolution leads to temporary, but significant increases in
potential energy, while we should only see a monotonic decrease along the free
deterministic flow. A more serious flaw is that EDGE does not appear to be a
saddle for the chosen quasipotential. For a saddle the potential should
decrease as we approach EDGE along its stable manifold on the basin boundary,
and further decrease when leaving along its unstable manifold. However, in our
case the dynamic enthalpy increases monotonically during the movement along
the basin boundary, as can be seen by the pseudo-trajectory along the basin
boundary during edge tracking (blue time series Fig. 11c). Hence, EDGE is
instead a local maximum of the dynamical enthalpy.

The M state we described in this section appears to be a robust steady state
of the system. Within statistical uncertainties, the same state is reached
when starting the edge tracking algorithm from two very close-by states on an
entirely different location of the basin boundary. In order to provide support
to this, we chose the intersection of the boundary and the straight line L1
(as opposed to L3), which is defined by a pair of different attractors
compared to the results presented above. Selected properties of the model
state obtained after 50 iterations of edge tracking using these initial
conditions are shown in Fig. 12. The shadowed states on the basin boundary
(red dotted line), as well as the estimated M state (black cross) match the
previous results very well (Fig. 12a), and the characteristic deep ocean
anomaly patterns are near-identical (compare Fig. 12b-h to Fig. 10a and Fig.
9).

Figure 12:  (a) Projection of the M states (green triangle and black cross)
and basin boundaries (red dotted and blue dashed lines) estimated by the edge
tracking procedure for two different initial conditions. The respective
initial conditions are taken by bisecting two very nearby trajectories (shown
in red and black) on opposite sides of the basin boundary shortly after they
separate. The two trajectories for the alternative initial conditions are
shown with reduced contrast. Also shown are the positions of the co-existing
attractors (colored symbols). (b) Anomaly of the depth-integrated dynamic
enthalpy of the M state with respect to the ON state. (c-h) Depth-integrated
temperature, salinity and density anomalies in the ocean below 1000 m of the M
state with respect to the ON and OFF attractors.

##  III Discussion

In this work we have investigated the global stability properties of the AMOC
by focusing, rather than on the competing attractive and competing circulation
regimes, on the properties of the dynamical saddles separating such states.

We harness the edge tracking algorithm to compute the M state, an unstable
AMOC configuration that lives on a basin boundary separating two attracting
regimes of a vigorous and collapsed circulation in a global ocean model. The
main physical relevance of such an unstable state is its potential role as
gateway of a noise-induced collapse of the AMOC. Determining its physical
characteristics uncovers the oceanographic changes to be expected leading up
to a noise-induced collapse. Since a priori the M state could be anywhere in
phase space, the AMOC strength could, for instance, first increase during an
ongoing transition to the OFF state, making an early detection non-obvious.
This specific scenario can be excluded with our results. Information on the M
state climatology can also be used to identify the physical mechanisms that
are necessary to drive a collapse, as well as specific perturbations that
would greatly increase the risk of a transition into the competing basin of
attraction. Furthermore, when the forcing parameter is increased towards the
tipping point where the reference attractor loses stability, fluctuations away
from the attractor are increasingly directed towards the M state (along the
instanton). This may inform us about observables, the fluctuations of which
could serve as physical (system-specific) early-warning signals.

By analyzing its climatology, we find that the M state lies consistently
between the ON and OFF states in terms of the upper ocean properties. Here it
features a freshening and cooling of the NA (and vice versa for the South
Atlantic), which is a known fingerprint of an AMOC slowdown. But the deep
ocean properties in EDGE are distinct and outside the range characterizing ON
and OFF. Most notably, EDGE features a cooling and freshening of almost the
entire Atlantic compared to both ON and OFF regimes. In oceanographic terms,
EDGE is less ’spicy’ (warm and salty) in the deep Atlantic compared to the
attractors. Since the OFF state lies between ON (most spicy) and EDGE (least
spicy), a transition from ON to EDGE and from ON to OFF both feature a
decreasing trend in spicyness. In contrast, in the deep ocean outside the
Atlantic there is a decrease in density when moving from ON to EDGE, which is
opposite to the trend when comparing ON and OFF. Thus, in terms of this
observable, a dangerous transition towards EDGE is the opposite as would be
expected when only knowing the ON-OFF anomaly and using the latter to
construct observables that signify an expected AMOC collapse. This shows the
importance of identifying all relevant characteristics of the M state (in
addition to OFF) in order to avoid a misinterpretation of oceanographic
changes leading up to an AMOC collapse. This is especially relevant as direct
observations of the AMOC strength have only been performed for the last two
decades, which is too short to detect changes considering the large
variability.

It may, however, be questioned whether deep ocean fingerprints are useful in
practice, due to difficulties in obtaining observations, and more
fundamentally because processes on very long time scales are likely involved
to accomplish the substantial change of the deep ocean properties (freshening
and cooling of the deep Atlantic beyond the values in OFF) required to reach
EDGE. It may thus be unlikely that a spontaneous collapse of the AMOC would
indeed occur as a passage via the instanton in the low-noise limit. Instead,
it may be a more likely scenario that, on shorter time scales, finite-
amplitude noise can kick the AMOC over the basin boundary in places within a
relatively short distance of the latter but far away from the M state. This
could be especially true if the noise has a considerable jump process
component (Lucarini  _et al._ , 2022). Due to the specific deep ocean
fingerprint, the likelihood of a transition via EDGE may also be quite
different for an AMOC collapse compared to an AMOC resurgence. If the AMOC
collapse is initiated by the cessation of NA deep water formation, as commonly
assumed, the Atlantic deep ocean properties would not immediately change
towards the EDGE conditions. This would need to happen gradually over time,
possibly by processes outside the NA. The AMOC resurgence, on the other hand,
may occur by a spontaneous onset of NA deep convection, which could in
principle flush the deep ocean with less spicy water much quicker and change
its properties towards EDGE.

The characteristics of the AMOC M state presented here need to be tested in
models with higher spatial resolution. Our model is a coarse-resolution ocean
model and does not resolve turbulent eddy dynamics, which may have an impact
on the mean ocean state and its stability. This could be tested using the same
ocean model in a higher resolution setup, but would require significantly
higher computational resources, both in order to perform the edge tracking and
to establish the stability landscape as an initial step. Further, since our
ocean-only model lacks atmospheric feedbacks, the coupled climate system could
in principle feature qualitatively different characteristics of the M state
and AMOC transition pathways, which possibly evade the long time scales
associated with the deep ocean. This needs to be tested with a coupled
atmosphere-ocean model, which is computationally expensive. It also remains to
be understood better whether the coupled climate system in the context of the
AMOC can indeed be considered as a system where the deterministic ocean is
forced by an effectively stochastic atmosphere that operates in the low-noise
limit such that the large-deviation theory applies, and as a result, where a
transition passes via the M state.

The Graham field theory indicates that an M state is a saddle in the
quasipotential landscape. While we have not derived a quantity that can serve
as quasipotential or Lyapunov function in the strict sense, EDGE indeed
features higher available potential energy compared to the attractors, as
measured by the relative magnitudes of the dynamic enthalpy. In a transition
between the attractors via EDGE, potential energy needs to be gained, which,
in the case of the M state we identified, needs to be achieved by bringing
relatively fresh (and thus light) water to the deep Atlantic. It would be
interesting to compare this specific transition path, obtained in a
deterministic way here, to actual noise-induced transitions paths when adding
atmospheric noise with finite strength to the system. Adding noise will also
shed light on whether the physical mechanism implied by the instanton implies
an asymmetric likelihood of ON$\rightarrow$OFF and OFF$\rightarrow$ON
transitions driven by short-term buoyancy changes. Further, empirical
quasipotentials obtained under stochastic forcing (Lucarini and Bódai, 2019,
2020; Margazoglou  _et al._ , 2021) may be compared to the measure of
potential energy employed here. By performing experiments at different values
of the control parameter, it could then be investigated how changes in EDGE,
and specifically the potential energy difference of attractors and EDGE, may
inform us of the likelihoods of spontaneous transitions, as suggested by
Arrhenius’ law. Changes in the edge state may also concern its dynamical
features, as opposed to the mean state discussed here, which could lead to
more useful fingerprints of an approaching AMOC collapse.

An important fundamental aspect that we have been unable to address here is
the study of the multitude of M states that may exist as our system features
several (at least seven) competing states for the chosen value of the control
parameter $F$. We have been able to only very partially explore the complexity
of the basin boundaries, each potentially supporting one or more M states,
which could in principle lead to a variety of distinct pathways for AMOC
tipping. The possibility of so-called Wada basins (Kennedy and Yorke, 1991),
where more than two basins of attraction share one boundary, makes an
investigation into the unknown number of co-existing M states even more
intriguing. We have so far not been able to find other M states that connect
the vigorous and collapsed AMOC regimes, since the only basin boundary we
could locate was between the two specific attractors ’ON’ and ’OFF’. But
additional M states may nevertheless exist and this should be investigated
with more extensive ensemble simulations or other techniques.

Acknowledgements We thank R. Nuterman and the Danish Center for Climate
Computing for supporting the simulations with the Veros ocean model.

Funding JL has received support from Danmarks Frie Forskningsfond under grant
no. 2032-00346B. VL has received support from the EU Horizon Europe project
ClimTIP (Grant No. 101137601), from the Marie Curie Action CRITICALEARTH
(Grant No. 956170), and from the EPSRC project LINK (Grant No. EP/Y026675/1).

Declaration of interests The authors report no conflict of interest.

Data availability statement The raw model simulation data that underlie the
findings of this study are deposited in the Electronic Research Data Archive
repository of the University of Copenhagen and can be accessed at
https://sid.erda.dk/sharelink/BMov9nurth.

Author ORCIDs J. Lohmann: https://orcid.org/0000-0002-6323-6243; V. Lucarini:
https://orcid.org/0000-0001-9392-1471.

Author contributions JL and VL conceived the study and designed the numerical
experiments. JL performed the numerical simulations and data analysis. JL and
VL wrote the paper.

## References

  * Westerhold  _et al._ (2020) T. Westerhold, N. Marwan, A. J. Drury, D. Liebrand, C. Agnini, E. Anagnostou, J. S. K. Barnet, S. M. Bohaty, D. D. Vleeschouwer, F. Florindo, T. Frederichs, D. A. Hodell, A. E. Holbourn, D. Kroon, V. Lauretano, K. Littler, L. J. Lourens, M. Lyle, H. Pälike, U. Röhl, J. Tian, R. H. Wilkens, P. A. Wilson, and J. C. Zachos, Science 369, 1383 (2020), https://www.science.org/doi/pdf/10.1126/science.aba6853 .
  * Ghil and Lucarini (2020) M. Ghil and V. Lucarini, Rev. Mod. Phys. 92, 035002 (2020).
  * Boers  _et al._ (2022) N. Boers, M. Ghil, and T. F. Stocker, Environ. Res. Lett. 17, 093006 (2022).
  * Rousseau  _et al._ (2023) D.-D. Rousseau, W. Bagniewski, and V. Lucarini, Scientific Reports 13, 11290 (2023).
  * Lenton (2013) T. M. Lenton, Annual Review of Environment and Resources 38, 1 (2013).
  * Alsulami and Petrovskii (2023) A. Alsulami and S. Petrovskii, Chaos, Solitons and Fractals 175, 114018 (2023).
  * Algeo and Shen (2024) T. J. Algeo and J. Shen, Natl Sci Rev 11, nwad237 (2024), 2053-714x Algeo, Thomas J Orcid: 0000-0002-3333-7035 Shen, Jun Orcid: 0000-0003-3759-6533 Journal Article Review China 2023/12/20 Natl Sci Rev. 2023 Sep 8;11(1):nwad237. doi: 10.1093/nsr/nwad237. eCollection 2024 Jan.
  * Lenton  _et al._ (2008) T. M. Lenton, H. Held, E. Kriegler, J. W. Hall, W. Lucht, S. Rahmstorf, and H. J. Schellnhuber, PNAS 105, 1786 (2008).
  * Armstrong McKay  _et al._ (2022) D. I. Armstrong McKay, A. Staal, J. F. Abrams, R. Winkelmann, B. Sakschewski, S. Loriani, I. Fetzer, S. E. Cornell, J. Rockström, and T. M. Lenton, Science 377, eabn7950 (2022).
  * Lucarini and Chekroun (2023) V. Lucarini and M. D. Chekroun, Nat Rev Phys 5, 744 (2023).
  * von der Heydt  _et al._ (2011) A. S. von der Heydt, P. Ashwin, C. D. Camp, M. Crucifix, H. A. Dijkstra, P. Ditlevsen, and T. M. Lenton, Glob. Planet. Change 197, 103399 (2011).
  * Imkeller and von Storch (2001) P. Imkeller and J. S. von Storch,  _Stochastic Climate Models_ (Birkhauser, Basel, 2001).
  * Hasselmann (1976) K. Hasselmann, Tellus XXVIII, 474 (1976).
  * Berner  _et al._ (2017) J. Berner, U. Achatz, L. Batté, L. Bengtsson, A. d. l. Cámara, H. M. Christensen, M. Colangeli, D. R. B. Coleman, D. Crommelin, S. I. Dolaptchiev, C. L. E. Franzke, P. Friederichs, P. Imkeller, H. Järvinen, S. Juricke, V. Kitsios, F. Lott, V. Lucarini, S. Mahajan, T. N. Palmer, C. Penland, M. Sakradzija, J.-S. von Storch, A. Weisheimer, M. Weniger, P. D. Williams, and J.-I. Yano, Bulletin of the American Meteorological Society 98, 565 (2017).
  * Saltzman (2001) B. Saltzman,  _Dynamical Paleoclimatology: Generalized Theory of Global Climate Change_ (Academic Press New York, New York, 2001).
  * Ashwin  _et al._ (2012) P. Ashwin, S. Wieczorek, R. Vitolo, and P. Cox, Phil. Trans. R. Soc. A 370, 1166 (2012).
  * Stommel (1961) H. Stommel, Tellus 13, 224 (1961).
  * Kuhlbrodt  _et al._ (2007) T. Kuhlbrodt, A. Griesel, M. Montoya, A. Levermann, M. Hofmann, and S. Rahmstorf, Rev. Geophys. 45, RG2001 (2007).
  * Weijer  _et al._ (2019) W. Weijer, W. Cheng, S. S. Drijfhout, A. V. Fedorov, A. Hu, L. C. Jackson, W. Liu, E. L. McDonagh, J. V. Mecking, and J. Zhang, Journal of Geophysical Research: Oceans 124, 5336 (2019), https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2019JC015083 .
  * Rahmstorf (2002) S. Rahmstorf, Nature 419, 207 (2002).
  * Liu  _et al._ (2020) W. Liu, A. V. Fedorov, S.-P. Xie, and S. Hu, Science Advances 6, eaaz4876 (2020), https://www.science.org/doi/pdf/10.1126/sciadv.aaz4876 .
  * Broecker  _et al._ (1990) W. S. Broecker, G. Bond, M. Klas, G. Bonani, and W. Wolfli, Paleoceanography 5, 469 (1990), https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/PA005i004p00469 .
  * Barker  _et al._ (2021) S. Barker, X. Zhang, L. Jonkers, S. Lordsmith, S. Conn, and G. Knorr, Paleoceanography and Paleoclimatology 36, e2020PA004200 (2021), e2020PA004200 2020PA004200, https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2020PA004200 .
  * Rahmstorf  _et al._ (2005) S. Rahmstorf, M. Crucifix, A. Ganopolski, H. Goosse, I. Kamenkovich, R. Knutti, G. Lohmann, R. Marsh, L. A. Mysak, Z. Wang, and A. J. Weaver, Geophys. Res. Lett. 32, L23605 (2005).
  * Hawkins  _et al._ (2011) E. Hawkins, R. S. Smith, L. C. Allison, J. M. Gregory, T. J. Woollings, H. Pohlmann, and B. de Cuevas, Geophys. Res. Lett. 38, L10605 (2011).
  * Boers (2021) N. Boers, Nat. Clim. Chang. 11, 680 (2021).
  * Ditlevsen and Ditlevsen (2023) P. Ditlevsen and S. Ditlevsen, Nat. Commun. 14, 4254 (2023).
  * van Westen  _et al._ (2024) R. M. van Westen, M. Kliphuis, and H. A. Dijkstra, Science Advances 10, eadk1189 (2024), https://www.science.org/doi/pdf/10.1126/sciadv.adk1189 .
  * Lohmann  _et al._ (2024) J. Lohmann, H. A. Dijkstra, M. Jochum, V. Lucarini, and P. D. Ditlevsen, Science Advances 10, eadi4253 (2024), https://www.science.org/doi/pdf/10.1126/sciadv.adi4253 .
  * Schneider  _et al._ (2008) T. M. Schneider, J. F. Gibson, M. Lagha, F. De Lillo, and B. Eckhardt, Phys. Rev. E 78, 037301 (2008).
  * Schneider and Eckhardt (2009) T. M. Schneider and B. Eckhardt, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 367, 577 (2009), https://royalsocietypublishing.org/doi/pdf/10.1098/rsta.2008.0216 .
  * Hof  _et al._ (2006) B. Hof, J. Westerweel, T. M. Schneider, and B. Eckhardt, Nature 443, 59 (2006).
  * Hof  _et al._ (2008) B. Hof, A. de Lozar, D. J. Kuik, and J. Westerweel, Phys. Rev. Lett. 101, 214501 (2008).
  * Skufca  _et al._ (2006) J. D. Skufca, J. A. Yorke, and B. Eckhardt, Phys. Rev. Lett. 96, 174101 (2006).
  * Vollmer  _et al._ (2009) J. Vollmer, T. M. Schneider, and B. Eckhardt, New Journal of Physics 11, 013040 (2009).
  * Lucarini and Bódai (2017) V. Lucarini and T. Bódai, Nonlinearity 30, R32 (2017).
  * Bódai and Lucarini (2020) T. Bódai and V. Lucarini, Chaos 30, 103105 (2020).
  * Lorenz (1975) E. Lorenz, GARP Publications series , 132 (1975).
  * Mehling  _et al._ (2024) O. Mehling, R. Börner, and V. Lucarini, Physica D: Nonlinear Phenomena 459, 134043 (2024).
  * Wieczorek  _et al._ (2011) S. Wieczorek, P. Ashwin, C. M. Luke, and P. M. Cox, Proc. R. Soc. A 467, 1243 (2011).
  * Lohmann  _et al._ (2021) J. Lohmann, D. Castellana, P. D. Ditlevsen, and H. A. Dijkstra, Earth Syst. Dynam. 12, 819 (2021).
  * Lohmann and Ditlevsen (2021) J. Lohmann and P. D. Ditlevsen, PNAS 118, e2017989118 (2021).
  * Graham (1987) R. Graham, in  _Fluctuations and Stochastic Phenomena in Condensed Matter_ (Springer, 1987) pp. 1–34.
  * Graham  _et al._ (1991) R. Graham, A. Hamm, and T. Tél, Phys. Rev. Lett. 66, 3089 (1991).
  * Touchette (2009) H. Touchette, Physics Reports 478, 1 (2009).
  * Schmeits and Dijkstra (2001) M. J. Schmeits and H. A. Dijkstra, J. Phys. Ocean. 31, 3435 (2001).
  * Plaut and Vautard (1993) G. Plaut and R. Vautard, J. Atmosph. Sci. 51, 210 (1993).
  * Glatzmaier and Roberts (1995) G. A. Glatzmaier and P. H. Roberts, Nature 377, 203 (1995).
  * Valet (2003) J.-P. Valet, Rev. Geophys. 41, 1004 (2003).
  * Bouchet  _et al._ (2019) F. Bouchet, J. Rolland, and E. Simonnet, Phys. Rev. Lett. 122, 074502 (2019).
  * Lucarini and Bódai (2019) V. Lucarini and T. Bódai, Phys. Rev. Lett. 122, 158701 (2019).
  * Lucarini and Bódai (2020) V. Lucarini and T. Bódai, Nonlinearity 33, R59 (2020).
  * Margazoglou  _et al._ (2021) G. Margazoglou, T. Grafke, A. Laio, and V. Lucarini, Proc. R. Soc. A 477, 20210019 (2021).
  * Börner  _et al._ (2024) R. Börner, O. Mehling, J. von Hardenberg, and V. Lucarini, Ocean Sciences Meeting 2024 , CC14C (2024).
  * Cessi (1994) P. Cessi, Phys. Oceanography 24, 1911 (1994).
  * Drijfhout  _et al._ (2013) S. Drijfhout, E. Gleeson, H. A. Dijkstra, and V. Livina, PNAS 110, 19713 (2013).
  * Kleppin  _et al._ (2015) H. Kleppin, M. Jochum, B. Otto-Bliesner, C. A. Shields, and S. Yeager, J. Climate 28, 7741 (2015).
  * Castellana and A. (2020) D. Castellana and D. H. A., Sci. Rep. 10, 20040 (2020).
  * Cini  _et al._ (2024) M. Cini, G. Zappa, F. Ragone, and S. Corti, npj Clim. Atm. Sci. 7, 31 (2024).
  * Häfner  _et al._ (2018) D. Häfner, R. L. Jacobsen, C. Eden, M. R. B. Kristensen, M. Jochum, R. Nuterman, and B. Vinter, Geosci. Model Dev. 11, 3299 (2018).
  * Häfner  _et al._ (2021) D. Häfner, R. Nuterman, and M. Jochum, JAMES 13, e2021MS002717 (2021).
  * Roberts  _et al._ (2013) C. D. Roberts, F. K. Garry, and L. C. Jackson, J. Clim. 26, 9155 (2013).
  * Rahmstorf  _et al._ (2015) S. Rahmstorf, J. E. Box, G. Feulner, M. E. Mann, A. Robinson, S. Rutherford, and E. J. Schaffernicht, Nature Clim Change 5, 475 (2015).
  * Caesar  _et al._ (2018) L. Caesar, S. Rahmstorf, A. Robinson, G. Feulner, and V. Saba, Nature 556, 191 (2018).
  * Jackson and Wood (2020) L. C. Jackson and R. A. Wood, J. Clim. 33, 7027 (2020).
  * Young (2010) W. R. Young, J. Phys. Ocean. 40, 394 (2010).
  * Nycander (2010) J. Nycander, J. Phys. Ocean. 41, 28 (2010).
  * Grafke and Vanden-Eijnden (2019) T. Grafke and E. Vanden-Eijnden, Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 063118 (2019), https://pubs.aip.org/aip/cha/article-pdf/doi/10.1063/1.5084025/13536944/063118_1_online.pdf .
  * Lucarini  _et al._ (2022) V. Lucarini, L. Serdukova, and G. Margazoglou, Nonlinear Process. Geophys. 29, 183 (2022).
  * Kennedy and Yorke (1991) J. Kennedy and J. A. Yorke, Physica D 51, 213 (1991).
  * Eden (2014) C. Eden, , pyOM2.0 Documentation, available at: https://wiki.cen.uni (2014).
  * Redi (1982) M. H. Redi, J. Phys. Oceanogr. 12, 1154 (1982).
  * Gent and McWilliams (1990) P. R. Gent and J. C. McWilliams, J. Phys. Oceanogr. 20, 150 (1990).
  * Gaspar  _et al._ (1990) P. Gaspar, Y. Grégoris, and J. M. Lefevre, J. Geophys. Res. 95, 16179 (1990).
  * Barnier (1998) B. Barnier, “Forcing the oceans,” in  _Ocean Modeling and Parameterization (ed. E. P. Chassignet and J. Verron)_ (Springer, 1998) pp. 45–80.
  * Uppala  _et al._ (2005) S. M. Uppala, P. W. Kållberg, A. J. Simmons, U. Andrae, V. D. C. Bechtold, M. Fiorino, J. K. Gibson, J. Haseler, A. Hernandez, G. A. Kelly, X. Li, K. Onogi, S. Saarinen, N. Sokka, R. P. Allan, E. Andersson, K. Arpe, M. A. Balmaseda, A. Beljaars, L. Berg, J. Bidlot, N. Bormann, S. Caires, F. Chevallier, A. Dethof, M. Dragosavac, M. Fisher, M. Fuentes, S. Hagemann, E. Hólm, B. J. Hoskins, L. Isaksen, P. Janssen, R. Jenne, A. P. Mcnally, J.-F. Mahfouf, J.-J. Morcrette, N. A. Rayner, R. W. Saunders, P. Simon, A. Sterl, K. E. Trenberth, A. Untch, D. Vasiljevic, P. Viterbo, and J. Woollen, Q. J. Roy. Meteor. Soc. 131, 2961 (2005).
  * Barnier  _et al._ (1995) B. Barnier, L. Siefridt, and P. Marchesiello, J. Marine Syst. 6, 363 (1995).
  * Amante and Eakins (2009) C. Amante and B. W. Eakins,  _US Department of Commerce, National Oceanic and Atmospheric Administration, National Environmental Satellite, Data, and Information Service, National Geophysical Data Center, Marine Geology and Geophysics Division Colorado_ , Tech. Rep. (2009).
  * Lohmann and Svensson (2022) J. Lohmann and A. Svensson, Clim. Past. 18, 2021 (2022).

##  Appendix A Details on the ocean model

We employ the primitive equation finite-difference global ocean model Veros, a
direct translation of the Fortran backend of the ocean model PyOM2 (Eden,
2014) into Python/JAX (Häfner  _et al._ , 2018, 2021). Mesoscale turbulence is
represented using the Redi (1982) and Gent and McWilliams (1990)
parameterization for isopycnal and thickness diffusion (diffusivity of 1 000
m2/s). We use the second order turbulence closure of Gaspar  _et al._ (1990)
to account for diapycnal mixing. Here, a background diffusivity of 10-5 m2/s
is employed, and static instabilities are removed by a gradual increase of the
vertical diffusivity as the instability is approached. This corresponds to a
smoothed variant of convective adjustment.

The heat exchange boundary condition is expressed by a first-order Taylor
expansion of the heat flux as a function of the anomaly of the modeled surface
temperature with respect to a fixed surface temperature climatology (Barnier,
1998). For this, we use ERA-40 climatologies (Uppala  _et al._ , 2005) of
surface temperature and heat flux, as well as a climatology of the derivative
of the heat flux with respect to changes in surface ocean temperature. The
latter is derived from ERA-40 following Barnier  _et al._ (1995) and has a
global yearly average 30.26 $WK^{-1}m^{-2}$. Wind stress forcing is also taken
from the ERA-40 reanalysis data. Freshwater exchanges with the atmosphere are
modeled by boundary conditions under which the sea surface salinity is relaxed
towards a present-day ERA-40 climatological field within a given relaxation
timescale. By choosing a long timescale of 2 years, oceanic salinity anomalies
are less efficiently damped by the atmospheric forcing compared to temperature
anomalies. This enables the positive salt advection feedback, which can lead
to AMOC multi-stability and tipping.

The bathymetry is obtained by smoothing the ETOPO1 global relief model (Amante
and Eakins, 2009) with a Gaussian filter to match the grid resolution. The
model domain ends at 80oN and thereby does not feature an Arctic connection of
ocean basins. The horizontal grid contains 90 longitudinal and 40 latitudinal
cells, where the latitudinal resolution increases from 5.3o at the poles to
2.1o at the Equator, as well as 40 vertical layers, which increase in
thickness from 23 m at the surface to 274 m at the bottom. For further model
details, see our previous studies that investigated AMOC tipping in this ocean
model (Lohmann and Ditlevsen, 2021; Lohmann and Svensson, 2022; Lohmann  _et
al._ , 2024).

In addition to the ERA-40 forcings, a freshwater anomaly is introduced in the
North Atlantic, in order to drive the system closer to the tipping point of a
collapsed AMOC. To this end, a salinity flux anomaly $\tilde{\phi}$ is applied
at the surface in the grid cells between 296oW to 0oW and 50oN to 75oN,
corresponding to an area of about $A=1.5\,\text{mio.}\,km^{2}$. Here we use
the equivalent total freshwater forcing $F=\tilde{\phi}AS_{ref}^{-1}$ as
control parameter, with the reference salinity $S_{ref}=3.5g\cdot kg^{-1}$.