# The Role of Sea-ice Processes on the Probability of AMOC Transitions

###### Abstract

Recent simulations performed with the Community Earth System Model (CESM) have
suggested a crucial role of sea-ice processes in AMOC hysteresis behaviour
under varying surface freshwater forcing. Here, we further investigate this
issue using additional CESM simulations and a novel conceptual ocean-sea-ice
box model. The CESM simulations show that the presence of sea ice gives rise
to the existence of statistical equilibrium states with a weak AMOC strength.
This is confirmed in the conceptual model, which captures the same AMOC
hysteresis behaviour as in the CESM simulation and where steady states are
computed versus forcing parameters. In the conceptual model, transition
probabilities between the different equilibrium states are determined using
rare event techniques. The transition probabilities from a strong AMOC state
to a weak AMOC state increase when considering sea-ice processes and indicate
that sea ice promotes these transitions. On the other hand, sea ice strongly
reduces the probabilities of the reverse transition from a weak AMOC state to
a strong AMOC state and this implies that sea ice also limits AMOC recovery.
The results here indicate that sea-ice processes play a dominant role in AMOC
hysteresis width and influence transition probabilities between the different
equilibrium states.

We develop a novel conceptual ocean-sea-ice box model to explain AMOC
hysteresis behaviour recently found in the Community Earth System Model and
determine how sea-ice processes influence the hysteresis width and the
probability of AMOC transitions.

##  1 Introduction

The Atlantic Meridional Overturning Circulation (AMOC) plays an important role
in regulating the global climate by redistributing heat and salinity through
the global ocean. The effects of anthropogenic climate change on the AMOC
provide reasons for concern as the AMOC is a tipping element in the climate
system (Armstrong McKay et al., 2022). There is an ongoing debate whether the
present-day AMOC is close to its tipping point (Boers, 2021; Ditlevsen and
Ditlevsen, 2023; Mehling et al., 2023), but it is known that AMOC tipping has
severe global climate impacts (Orihuela-Pinto et al., 2022; van Westen et al.,
2024). So far, historical AMOC reconstructions indicate that its strength is
declining since 1900 (Caesar et al., 2018, 2021) and projections of climate
models indicate that this decline will continue under future climate change
and its decline is insensitive to the used emission scenario (Weijer et al.,
2020).

The potential tipping of the AMOC has been studied in a hierarchy of models
(Dijkstra, 2023), from simple box models (Stommel, 1961; Cessi, 1994;
Cimatoribus et al., 2014), ocean-only models (Dijkstra, 2007), to various
Earth system Models of Intermediate Complexity (EMICs) (Rahmstorf et al.,
2005), early Global Climate Models (GCMs) (Hawkins et al., 2011; Hu et al.,
2012), and in modern Earth System Models (ESMs) (Jackson et al., 2022). In box
models and ocean-only models, steady state calculations have shown that multi-
stable AMOC regimes can exist. Transitions between strong and collapsed AMOC
states are caused by the salt-advection feedback, where a freshwater anomaly
in the North Atlantic weakens the AMOC strength and reduces the associated
salinity transport. This response then amplifies the original freshwater
anomaly (Marotzke, 2000). In EMICs and ESMs, an impression of such a multi-
stable regime can be obtained by slowly increasing the surface freshwater
forcing to induce an AMOC collapse. When reversing the surface freshwater
forcing, the AMOC often recovers at smaller values of this forcing compared to
the collapse, giving rise to hysteresis behaviour (Rahmstorf et al., 2005;
Hawkins et al., 2011; Hu et al., 2012).

The presence of a multi-stable AMOC in a hierarchy of models, and possibly in
the actual climate system, has far reaching consequences. When the present-day
AMOC is in such a regime, there is a finite probability of a noise-induced
transition between the different equilibrium states, given the presence of
stochastic surface forcing. This probability will increase when approaching
the tipping point; in addition, rate-induced transitions may occur Ashwin et
al. (2012); Lohmann and Ditlevsen (2021). Transition probabilities between the
different AMOC states, which may arise under stochastic surface freshwater
flux noise, have so far been only studied in idealised models (Castellana et
al., 2019; Baars et al., 2021). For state-of-the-art climate models
determining transition probabilities is problematic since this is
computationally very expensive, in particular when the transition
probabilities are very low (Jacques-Dumas et al., 2023).

This study is motivated by recent results where such hysteresis behaviour was
found in the Coupled Model Intercomparison Project phase 5 (CMIP5) version of
the Community Earth System Model (CESM) (van Westen and Dijkstra, 2023; van
Westen et al., 2024). The transition (‘collapse’) between the strong (‘AMOC
on’) to the weak (‘AMOC off’) state was very similar to the one expected from
ocean-only models (Dijkstra, 2007), whereas the reverse transition
(‘recovery’) behaved quite differently. It was shown that sea-ice ocean
interactions play an important role in the AMOC recovery because the
relatively weak ($<$5 Sv, 1 Sv = 106 m3 s-1) northward overturning cell in the
Atlantic Ocean is strongly suppressed by the extensive North Atlantic sea-ice
cover (van Westen and Dijkstra, 2023). However, the precise effects of the sea
ice on the multi-stable AMOC regime could not be investigated in CESM, where
the surface freshwater forcing was transient.

In this paper, we show that the sea ice actually increases the hysteresis
width in the CESM and explain the physics of this result using a novel
conceptual ocean-sea-ice model. Using this conceptual model, we also extend
the study of Castellana et al. (2019) by looking at transition probabilities
between the different equilibrium AMOC states in the presence of sea ice. The
conceptual ocean-sea-ice model is described in Section 2. Next in Section 3,
we first present additional CESM results and then analyse the steady states
and transitions between the steady states in the conceptual model. Finally in
Section 4, we summarise and discuss the results.

##  2 Models and Methods

###  2.1 CESM configuration

In van Westen and Dijkstra (2023), version 1 of the CESM was used, with
horizontal resolutions of 1∘ for the ocean/sea ice and 2∘ for the
atmosphere/land components. The CESM simulation was initialised from the end
of a 2,800-year long pre-industrial control simulation (Baatsen et al., 2020)
under constant greenhouse gas, solar and aerosol forcings set to pre-
industrial levels. The hysteresis behaviour in the CESM (van Westen and
Dijkstra, 2023) was obtained by varying the surface freshwater forcing
($F_{H}$) between latitudes 20∘N and 50∘N in the Atlantic Ocean. This
freshwater flux anomaly was compensated over the rest of the domain to
conserve salinity, this a typical set-up for a hosing experiment. The surface
freshwater forcing linearly increased in time with a rate of $3\times 10^{-4}$
Sv yr-1 up to a maximum of $F_{H}=0.66$ Sv in model year 2,200. This rate of
change in $F_{H}$ is quite small with the aim to stay close to equilibrium,
often referred to as a so-called quasi-equilibrium approach (Rahmstorf et al.,
2005; Hu et al., 2012). At the end of the simulation the forcing was reversed
and $F_{H}$ decreased at the same rate and reaches zero in model year 4,400.

The different AMOC states in van Westen and Dijkstra (2023) are not in
statistical equilibrium (Hawkins et al., 2011), given the quasi-equilibrium
approach, and therefore we performed additional CESM simulations under fixed
freshwater forcing values. From the hysteresis simulation we branch off 4
simulations from model years 600 and 3,800 ($F_{H}=0.18$ Sv) and from model
years 1,500 and 2,900 ($F_{H}=0.45$ Sv), and keep the value of $F_{H}$
constant during each simulation. The AMOC collapse (recovery) occurs around
model year 1,750 (4,100), this motivates the branching points from
$F_{H}=0.45$ Sv) ($F_{H}=0.18$ Sv). The branched simulations are integrated
for 500 years and we determine their statistical properties over the last 50
years.

###  2.2 A conceptual ocean-sea-ice model

The new model is based on the ocean box model which was developed by
Cimatoribus et al. (2014) and extended by Castellana et al. (2019) (hereafter
called CCM, see Figure 1). The CCM consists of conservation equations for the
salinity content of boxes s, ts, t and n, one equation to conserve salinity
(from which the salinity of box d is derived) and one equation for the time-
varying pycnocline depth ($D$). The AMOC strength in the northern box
($q_{N}$) is given by

$q_{N}=\eta\frac{\rho_{n}-\rho_{ts}}{\rho_{0}}D^{2},$ (1)

where $\eta$ is a downwelling (pumping) constant, $\rho_{n}-\rho_{ts}$ is the
density difference between box n and box ts and $\rho_{0}$ is a reference
density. The densities are determined from a linear equation of state. In the
CCM, the temperatures for box n and box ts are input parameters (5∘C and 10∘C,
respectively) and are fixed in time. The main forcing parameter is the
asymmetric freshwater flux ($E_{A}$) from box s to box n. Note that the hosing
set-up in the CCM is slightly different than in the CESM where fresh water is
added between 20∘N and 50∘N and compensated elsewhere. Under large $E_{A}$,
the density difference $\rho_{n}-\rho_{ts}$ becomes negative and the
circulation reverses (red dashed arrows become red dotted arrows in Figure 1).
The CCM is extensively discussed in Castellana et al. (2019) and, as we use it
for comparison to the new model where sea ice is included, the governing
equations and input parameter settings are provided in the Appendix (see
equations 6).

Figure 1: Schematic representation of the 5-box AMOC model, adapted from
Cimatoribus et al. (2014) and Castellana et al. (2019) (the CCM). The red
arrows represent the volume fluxes, whereas the dashed and dotted arrows
indicate the AMOC on and AMOC off states, respectively. The cyan and blue
arrows are the gyre transport and freshwater fluxes. The brown arrows are the
heat fluxes with the overhead atmosphere for each surface box (i.e., box s,
ts, t and n).

Our aim is to couple sea-ice processes to the CCM and hence we first extend it
by including oceanic temperatures. The governing equations for temperature
(see equations 7 in the Appendix) are very similar to the ones for salinity.
The main difference is that the equations for the surface boxes (box s, ts, t
and n) have one additional term which represents the heat exchange with the
atmosphere (brown arrows in Figure 1). The change in ocean heat content for
each surface box ($i$ = s, ts, t and n) and box d is then given by (under the
assumption of constant heat capacity and reference density):

$\displaystyle\frac{\mathrm{d}(V_{i}T_{i})}{\mathrm{d}t}$
$\displaystyle=F_{i}-\lambda^{a}A_{i}(T_{i}-T_{i}^{a}),$ (2a)
$\displaystyle\frac{\mathrm{d}(V_{\mathrm{d}}T_{\mathrm{d}})}{\mathrm{d}t}$
$\displaystyle=F_{\mathrm{d}}.$ (2b)

Here, $F_{i}$ is the advective oceanic heat transport, $V_{i}$ the volume,
$A_{i}$ the horizontal surface area, $\lambda^{a}$ the heat exchange rate with
the atmosphere, $T_{i}$ the ocean temperature and $T_{i}^{a}$ the overhead
atmospheric temperature. Note that the temperature for box d is a dynamical
quantity whereas the salinity of this box is diagnostically derived from
salinity conservation.

The heat exchange rate $\lambda^{a}$ and atmospheric temperatures are
additional parameters in the extended CCM (hereafter E-CCM) that need to be
tuned (Appendix and Figure 9). For an arbitrary value of $E_{A}$ (and for the
AMOC on state), we tune these parameters such that the ocean temperatures for
box n and box ts are 5∘C and 10∘C, respectively. We mainly use
$\lambda^{a}=3.5\times 10^{-6}$ m s-1 below and for the given $E_{A}$ the AMOC
on state of the E-CCM is identical to that of the CCM. We used $E_{A}=0.0644$
Sv for tuning and varying this tuning value slightly shifts the bifurcation
points in the E-CCM. In the CCM this is the minimum value of $E_{A}$ for which
the freshwater transport carried by the AMOC at 30∘S, below indicated by
$F_{\mathrm{ovS}}$, becomes negative (Castellana et al., 2019). The quantity
$F_{\mathrm{ovS}}$ is determined from the model as

$F_{\mathrm{ovS}}=-\frac{q_{S}}{S_{0}}\left(S_{\mathrm{ts}}-S_{\mathrm{d}}\right),$
(3)

where $S_{0}$ is a reference salinity. The quantity $F_{\mathrm{ovS}}$ is
important as it is a measure for the salt-advection feedback strength. The
role of $F_{\mathrm{ovS}}$ is more versatile in the CCM because when it is
negative (positive) in the AMOC on state the model is in its multi-stable
(monostable) regime.

Finally, the sea-ice processes are represented as follows and specifically for
$\lambda^{a}=3.5\times 10^{-6}$ m s-1. When the temperature for box n is lower
than 5∘C (its reference temperature), we assume that the sea-ice fraction of
box n ($f_{\mathrm{n}}$) grows linearly by 20% per degree cooling (Figure
10a). This is an arbitrary choice for the sea-ice growth and we do not use the
freezing temperature for sea-ice growth as box n represents the volume-
averaged temperature over the higher latitudes. To represent now the effect
that sea ice suppress upper ocean mixing (Lin et al., 2023), we multiply the
AMOC strength by a given AMOC reduction factor $0\leq r_{q}\leq 1$ such that
the AMOC strength becomes

$q_{N}^{\mathrm{ice}}=r_{q}~{}q_{N},$ (4)

where $q_{N}$ is the AMOC strength from (1). Since a continuous step-like
change of the AMOC is found in the CESM (see Figure 4 in van Westen and
Dijkstra (2023)), we use a similar function for the reduction factor (Figure
10b). We choose

$r_{q}=\frac{1-r_{q}^{\mathrm{min}}}{2}\tanh(0.1\times(0.5-f_{\mathrm{n}}))+\frac{1+r_{q}^{\mathrm{min}}}{2}.$
(5)

such that $r_{q}\sim 1$ when the sea-ice fraction $f_{\mathrm{n}}<0.4$ (40%)
and, for $f_{\mathrm{n}}>0.6$ (60%), $r_{q}$ decreases to a minimum value
$r_{q}^{\mathrm{min}}$. When $r_{q}=1$ is chosen, the model does not represent
any sea ice effects.

In the CESM results (van Westen and Dijkstra, 2023), there is a factor 10
difference in spatially-averaged mixed layer depth over the Irminger basin
between the two AMOC states (which would suggest $r_{q}^{\mathrm{min}}\approx
0.1$) with local differences up to a factor 20 difference
($r_{q}^{\mathrm{min}}\approx 0.05$). We will vary $r_{q}^{\mathrm{min}}$
below to investigate its sensitivity in the the E-CCM. Note that the value of
$r_{q}$ depends on the sea-ice fraction and hence on the temperature
$T_{\mathrm{n}}$, but sea ice is not dynamically resolved in the E-CCM and its
effect is only parameterised.

###  2.3 Steady States and Transition Probabilities

To determine the bifurcation diagrams of the CCM and E-CCM we implemented the
model equations in the continuation software AUTO-07p (Doedel et al., 2007,
2021). The AUTO software is able to directly compute steady states versus
model parameters without time integration by using a pseudo-arclength
continuation combined with a Newton-Raphson method. An additional advantage of
AUTO is its ability to detect special points such as Hopf and saddle-node
bifurcations, which are both present in the CCM and E-CCM. The accuracy of
AUTO is determined by three coefficients representing the absolute and
relative accuracy of the solution, and the accuracy for locating special
points. All these coefficients are set on a minimum value of 10-6 and a higher
accuracy is used where necessary to determine an accurate bifurcation diagram.

In the multi-stable regime, we study noise-induced transitions in the
stochastic version of the E-CCM. Stochastic variations in the asymmetric
freshwater flux forcing are considered by adding a term
$f_{\sigma}E_{A}S_{0}\zeta(t)$ to the salinity equation of box n (and
subtracting it from the salinity equation of box s), where $\zeta(t)$ is a
zero mean, delta time correlated (white) noise process and $f_{\sigma}$
controls the strength of the noise. Following Castellana et al. (2019), we use
the Trajectory Adaptive Multilevel Splitting (TAMS) algorithm (Lestang et al.,
2018) to determine probabilities of transitions between different steady AMOC
states in the E-CCM within a given time interval $T$.

##  3 Results

###  3.1 Statistical Steady States in CESM

The branched simulations from the hysteresis experiment identify multiple
statistical steady AMOC states in the CESM (Figure 2) for both $F_{H}=0.18$ Sv
and $F_{H}=0.45$ Sv.

The first state is obtained from the strong northward overturning state for
$F_{H}=0.18$ Sv (branched from model year 600). The AMOC strength appears
statistically stationary after 500 years of the branched simulation. When we
consider the time means (450 – 500 years after the branching point), we find
an AMOC strength and $F_{\mathrm{ovS}}$ of 15.7 Sv and 0.07 Sv, respectively.
The results from the quasi-equilibrium simulation (model years 575 – 624) are
close to the latter equilibrium values with 14.8 Sv for AMOC strength and 0.10
Sv for $F_{\mathrm{ovS}}$. For $F_{H}=0.18$ Sv (branched from model year
3,800), there is a second statistical steady state. This state has a weak and
shallow ($<$1,000 m) northward overturning cell (van Westen and Dijkstra,
2023) in the Atlantic Ocean with an equilibrium AMOC strength of 4.2 Sv and
$F_{\mathrm{ovS}}$ of $-$0.10 Sv. The quasi-equilibrium values (model years
3,775 – 3,824) are 3.2 Sv and $-$0.07 Sv for the AMOC strength and
$F_{\mathrm{ovS}}$, respectively.

For $F_{H}=0.45$ Sv (branched from model year 1,500), the AMOC is still in the
northward overturning regime, but its strength has declined under the stronger
freshwater forcing. The equilibrium values are 12.4 Sv for AMOC strength and
$-$0.16 Sv for $F_{\mathrm{ovS}}$, while the quasi-equilibrium approach (model
years 1,475 – 1,524) values are 11.9 Sv and $-$0.10 Sv, respectively. As
discussed in van Westen et al. (2024), the sign and variance in
$F_{\mathrm{ovS}}$ are good indicators of the approach of the AMOC tipping
point. The $F_{\mathrm{ovS}}$ variance for $F_{H}=0.45$ Sv ($2.1\times
10^{-4}$ Sv2) is indeed larger than at $F_{H}=0.18$ Sv ($0.33\times 10^{-4}$
Sv2), which is directly visible from the branched time series (Figure 2d). A
linear trend was removed before determining the variances. The larger variance
indicates that the AMOC loses resilience and is closer to its tipping point.
The variances in AMOC strength are quite similar for both branches and are
about 0.3 Sv2. For $F_{H}=0.45$ Sv another statistical steady state exists;
this is the collapsed AMOC state (branched from model year 2,900). This
equilibrium state has an AMOC strength of 0.08 Sv and $F_{\mathrm{ovS}}$ of
0.11 Sv and the values for the quasi-equilibrium approach (model years 2,875 –
2,924) are close to the equilibrium values.

The quasi-equilibrium simulation remains reasonably close to the steady states
but also indicates that the freshwater forcing changes are faster than oceanic
adjustment time scales. Ideally one would like to change the freshwater
forcing even slower (than $3\times 10^{-4}$ Sv yr-1) and determine more steady
states, preferably closer to the tipping points, but these simulations require
even more computational resources and are at the moment not feasible. As
argued in van Westen and Dijkstra (2023), one would expect the AMOC recovery
to occur at the value of $F_{H}$ where $F_{\mathrm{ovS}}$ changes sign on the
AMOC on branch (black curve in Figure 2b), that is around year 3,400 in the
quasi-equilibrium approach. The results in Figure 2 show instead that a
statistical steady AMOC off state still exists for $F_{H}=0.18$ Sv and hence
the late recovery at year 4,100 is not due to an overshoot but to the fact
that, in a model with sea ice, AMOC off states can apparently extend to
positive values of $F_{\mathrm{ovS}}$. We will explore this below in more
detail using the conceptual model.

Figure 2: (a): The AMOC strength at 1,000 m and 26∘N. The cyan-coloured lines
in panels a and b indicate the magnitude of the freshwater forcing $F_{H}$ and
the yellow shading in all panels indicates observed ranges for the presented
quantity (Garzoli et al., 2013; Mecking et al., 2017; Smeed et al., 2018;
Worthington et al., 2021). The statistical equilibria (i.e., steady states in
blue, panels a and b) of model years 600, 1,500, 2,900 and 3,800 are also
shown, where the marker indicates the mean and the error bars show the minimum
and maximum over the last 50 years of the 600-year long branched simulations.
Inset: The hosing experiment where fresh water is added to the ocean surface
between 20∘N – 50∘N in the Atlantic Ocean ($+F_{H}$) and is compensated over
the remaining ocean surface ($-F_{H}$). The black sections indicate the 26∘N
and 34∘S latitudes over which the AMOC strength and $F_{\mathrm{ovS}}$ are
determined, respectively. (b): The freshwater transport by the AMOC at 34∘S,
$F_{\mathrm{ovS}}$. (c & d): Similar to panels a,b, but now the entire
branched simulations are displayed. The arrows indicate the beginning of the
branches and the direction of time for each branch.

###  3.2 Steady States in the E-CCM

To investigate the important role of sea ice in the AMOC hysteresis, we
compute steady states in the deterministic E-CCM. First we analyse the case
without sea ice (by setting $r_{q}^{\mathrm{min}}=1$) and vary the asymmetric
freshwater flux forcing ($E_{A}$). The bifurcation diagrams for this case are
shown in Figure 3a for two values of $\lambda^{a}$. Both bifurcation points
are saddle node bifurcations, similar to the CCM, and close to the saddle node
bifurcations ($\Delta E_{A}=0.002$ Sv) there are (subcritical) Hopf
bifurcations which give unstable periodic orbits (not shown). The position of
the bifurcation points is dependent on the heat exchange rate $\lambda^{a}$
with the atmosphere and lower values of $\lambda^{a}$ (less restoring) result
in a larger value of $E_{A}$ at the saddle-node bifurcation points. In the
limit of $\lambda^{a}\rightarrow\infty$, we find the same bifurcation points
as in the CCM (dashed lines in Figure 3d). The hysteresis width, that is the
distance between $E_{A}^{1}$ and $E_{A}^{2}$, is about 0.3 Sv and does not
depend much on $\lambda^{a}$.

To understand this behaviour of the shift of $E^{1}_{A}$ with $\lambda^{a}$,
we analyse the AMOC strength response under varying $E_{A}$. The AMOC strength
is represented by $q_{N}$ and is proportional to the density difference
between box n and box ts ($\Delta\rho$), and, given the linear equation of
state, it boils down to $\Delta\rho\sim\Delta T-\Delta S$. Under the
increasing freshwater forcing the salinity difference becomes larger and this
decreases the density difference resulting in AMOC weakening (Figure 3a). The
temperature of box n ($T_{\mathrm{n}}$, Figure 3b) decreases as this box
receives less heat from box t under AMOC weakening, which enhances the
temperature difference between box n and box ts (the temperature of box ts
remains fairly constant). This larger temperature difference between box n and
box ts increases the density difference, so there are two counteracting
effects that control the density difference. The combined effect of the
temperature and salinity responses under the freshwater forcing shifts the
saddle-node bifurcation point ($E_{A}^{1}$, Figure 3d) to higher values of
freshwater forcing when comparing this to the standard CCM (dashed lines in
Figure 3d). The stabilising temperature effect is larger at smaller
$\lambda^{a}$ (less restoring) and hence a larger value of $E^{1}_{A}$ is
needed to destabilise the AMOC (Figure 3) at the saddle-node bifurcation.

Figure 3: (a – c): The steady state solutions for the E-CCM (no sea ice) for
the (a): AMOC strength, (b): temperature of box n and (c): $F_{\mathrm{ovS}}$.
The red and blue curves indicate the steady AMOC on and AMOC off state,
respectively, and for two different values for $\lambda^{a}$ (solid and dashed
lines). The (dash)dotted curves are the unstable branches and the saddle-node
bifurcation points are indicated by the circled markers. (d): The bifurcation
points as a function of freshwater forcing for varying $\lambda^{a}$, the
dashed lines indicate the bifurcation points for the CCM.

When the AMOC collapses, both temperature and salinity (strongly) decline in
box n, these changes are much smaller for box ts. The density gradient is
influenced by both salinity and temperature, the temperature response partly
compensates the salinity response and hence the bifurcation point
($E_{A}^{2}$, Figure 3d) also shifts to higher values of $E_{A}$ compared to
that of the CCM. Note that the temperature of box n does not vary much under
the varying freshwater forcing when AMOC is in its off state (Figure 3b). This
box is in equilibrium as it only receives heat from its overhead atmosphere
and from the gyre transport with box t. The Hopf bifurcations near $E_{A}^{1}$
and $E_{A}^{2}$ have associated periods of about 60 to 200 years and 3,000 to
6,000 years, respectively. For relatively large values of $\lambda^{a}$, the
value of $E_{A}$ where $F_{\mathrm{ovS}}$ changes sign on the AMOC on branch
is close to $E_{A}^{2}$, similar as in the CCM. But for small values of
$\lambda^{a}$ (e.g., $\lambda^{a}=3.5\times 10^{-6}$ m s-1), this is clearly
not the case. In other words, negative $F_{\mathrm{ovS}}$ does not exclusively
indicate the multi-stable regime.

Next, the case with sea-ice representation in the E-CCM is considered for
$\lambda^{a}=3.5\times 10^{-6}$ m s-1. Figure 4 shows the bifurcation diagram
for $r_{q}^{\mathrm{min}}=0.02$, where the saddle node bifurcation $E^{2}_{A}$
in Figure 3 has moved to $E^{3}_{A}$. The saddle-node $E_{A}^{1}$ is hardly
affected as the AMOC reduction factor $r_{q}$ is about 1 (relatively low sea-
ice fractions when $f_{\mathrm{n}}<35\%$) when in the AMOC on state. At the
value $E_{A}^{2}$ there is no longer a saddle-node bifurcation, but it remains
an important point as $q_{N}$ switches sign here; in the case without sea ice
this would then result in AMOC recovery. However, in the presence of sea ice,
the AMOC reduction factor is very strong (here $r_{q}\approx 0.02$) and a new
branch of equilibrium states appears in which the AMOC remains weak (less than
3 Sv). For decreasing values of the freshwater forcing, the AMOC becomes
stronger and the sea-ice fraction starts declining until the AMOC reaches
about 3 Sv at the new bifurcation point $E_{A}^{3}$. There are again
(subcritical) Hopf bifurcations close to the saddle nodes $E_{A}^{1}$ and
$E_{A}^{3}$ with associated periods of 50 to 250 years and 200 to 1,250 years,
respectively.

Figure 4: (a – c): The steady state solutions for the E-CCM for the (a): AMOC
strength, (b): temperature of box n and sea-ice fractions and (c):
$F_{\mathrm{ovS}}$. The red, blue and black curves indicate the AMOC on, AMOC
off and AMOC weak state, respectively, for $\lambda^{a}=3.5\times 10^{-6}$ m
s-1 and $r_{q}^{\mathrm{min}}=0.02$. The (dash)dotted curves are the unstable
branches and the saddle-node bifurcation points are indicated by the circled
markers, the diamond marker is not a saddle-node bifurcation point. The inset
in panel c is a zoom in around the $F_{\mathrm{ovS}}$ minimum of the AMOC on
branch. (d): The bifurcation points as a function of $E_{A}$ for varying
$r_{q}^{\mathrm{min}}$ and for $\lambda^{a}=3.5\times 10^{-6}$ m s-1. The two
insets shows the steady state for AMOC strength ($q_{N}^{\mathrm{ice}}$,
similar to panel a) for $r_{q}^{\mathrm{min}}=0.038$ and
$r_{q}^{\mathrm{min}}=0.054$, centred around the weak state.

When decreasing the freshwater forcing from $E_{A}^{2}$, $F_{\mathrm{ovS}}$
becomes slightly positive and then changes sign. In this freshwater forcing
range, the AMOC strength is very weak ($q_{N}^{\mathrm{ice}}<0.6$ Sv) and
$q_{S}$ remains negative. Lowering the freshwater forcing further results in a
stronger AMOC and $q_{S}$ becomes positive, which induces a negative
$F_{\mathrm{ovS}}$. In the CCM, negative values of $F_{\mathrm{ovS}}$ in the
AMOC on state mark the multi-equilibrium regime (Castellana et al., 2019) and
here for $r_{q}^{\mathrm{min}}=0.02$ this is (by coincidence) also the case
(Figure 4c). However, for $r_{q}^{\mathrm{min}}<0.02$ the bifurcation point
$E_{A}^{3}$ shifts to lower values of the freshwater forcing and then positive
$F_{\mathrm{ovS}}$ are in multi-equilibrium regime, similar as in the CESM.
The quantity $F_{\mathrm{ovS}}$ is still an important indicator in the CESM as
its sign change coincides with the initial AMOC recovery but the AMOC strength
is strongly suppressed by sea-ice processes (van Westen and Dijkstra, 2023).
When using a larger $\lambda^{a}$ and a different sea-ice response function in
the E-CCM, we can tune the model such that $E_{A}^{2}$ is close the
$F_{\mathrm{ovS}}$ sign change, similar as in the CESM.

It is interesting that the weak state splits into two sub-branches for
$0.030\leq r_{q}^{\mathrm{min}}\leq 0.052$ (inset in Figure 4d for
$r_{q}^{\mathrm{min}}=0.038$); the AMOC on and AMOC off states remain very
similar. The two sub-branches differ in $F_{\mathrm{ovS}}$ sign and there are
subcritical Hopf bifurcations near three saddle-node bifurcations located on
$E_{A}^{1}$, $E_{A}^{2}$, and $E_{A}^{3}$. The associated periods for the Hopf
bifurcations near $E_{A}^{1}$ and $E_{A}^{3}$ are 50 – 500 years and 150 – 500
years, respectively. The Hopf bifurcation near $E_{A}^{2}$ has a much longer
associated period of 44,000 – 52,000 years. For
$0.052<r_{q}^{\mathrm{min}}\leq 0.060$ only the weak state with negative
$F_{\mathrm{ovS}}$ remains (inset in Figure 4d for
$r_{q}^{\mathrm{min}}=0.054$) and $E_{A}^{2}$ is now a saddle node (similar to
Figure 3). For increasing $r_{q}^{\mathrm{min}}$, the saddle node $E_{A}^{3}$
shifts to larger values of $E_{A}$ and has a maximum AMOC strength of 3 Sv.
The AMOC weak state only appears when the AMOC reduction factor is
sufficiently low ($r_{q}^{\mathrm{min}}\leq 0.060$) and for very low values of
this factor ($r_{q}^{\mathrm{min}}<0.015$) a negative freshwater forcing is
required to return to the AMOC on state (Figure 4d). Note that for
$r_{q}^{\mathrm{min}}>0.060$ the bifurcation point $E_{A}^{3}$ disappears and
then the results are similar to Figure 3. The hysteresis width also varies for
varying $r_{q}^{\mathrm{min}}$. For values of $r_{q}^{\mathrm{min}}\geq
0.046$, it is a constant width of 0.3 Sv while the hysteresis width increases
for lower values for $r_{q}^{\mathrm{min}}$. For example, for
$r_{q}^{\mathrm{min}}=0.02$ it is about 0.4 Sv which is comparable to that of
the CESM.

###  3.3 Quasi-equilibrium Experiment in the E-CCM

When varying the freshwater flux forcing over time with a rate of $3\times
10^{-4}$ Sv yr-1 (similar to the CESM), we find qualitatively the same
hysteresis behaviour in the E-CCM (Figure 5) as in the CESM (Figure 2). Note
that we continued the simulation beyond model year 4,400 under no freshwater
forcing to run the model into equilibrium. We find an AMOC tipping event
around model year 1,600 and $F_{\mathrm{ovS}}$ goes through a minimum and
further equilibrates to the reversed circulation. The steady states show the
$F_{\mathrm{ovS}}$ minimum prior to the AMOC tipping point (inset in Figure
4c) with a difference of only $\Delta E_{A}=8\times 10^{-4}$ Sv. Other
idealised models find (qualitatively) a similar timing between
$F_{\mathrm{ovS}}$ minimum and AMOC tipping (Dijkstra, 2007; Huisman et al.,
2010) and in the CESM (van Westen et al., 2024).

The AMOC remains in its off state until around model year 3,900 and
$q_{N}^{\mathrm{ice}}$ becomes positive and $F_{\mathrm{ovS}}$ switches sign.
Eventually the AMOC transitions to the on state and strongly overshoots (to 68
Sv). This overshoot is a model artefact as the AMOC strength is then
determined by the (instantly felt) density difference between box n and box ts
and by the pycnocline depth which is about twice as deep. The AMOC in the CESM
is not only controlled by this density difference and the pycnocline depth,
which explains why the CESM does not have such a large overshoot. There is a
second overshoot in AMOC strength of about 22.9 Sv (around model year 4,400),
which is more realistic and dynamically driven. Advection of temperature and
salinity controls the adjustment time scale from the weak AMOC to the AMOC on
state and this takes about a few centuries, which is much faster than in the
CESM (or real ocean), where vertical mixing also plays an important role.
Nevertheless, prior to the full adjustment to the AMOC on state, the different
temperature and salinity gradients (w.r.t. the initial state) allow for a
temporarily stronger AMOC strength. At the same time, $F_{\mathrm{ovS}}$
reaches its minimum and then switches sign when the salinity in the boxes
further adjust to the AMOC on state.

Figure 5: (a): The AMOC strength ($q_{N}^{\mathrm{ice}}$) and (b): the
freshwater transport by the AMOC at 30∘S in the E-CCM for
$\lambda^{a}=3.5\times 10^{-6}$ m s-1 and $r_{q}^{\mathrm{min}}=0.02$. The
cyan-coloured lines in panels a and b indicate the magnitude of $E_{A}$. The
yellow shading in all panels indicates observed ranges for AMOC strength and
$F_{\mathrm{ovS}}$. Note that after model year 4,400 the freshwater forcing is
held constant at $E_{A}=0$ Sv.

The transient responses in the E-CCM are comparable to that of the CESM. The
deviations in the quasi-equilibrium approach w.r.t. the steady states of the
AMOC on branch become larger when approaching the tipping point. Note that the
timing of AMOC collapse and recovery is not identical to that of the CESM, it
is possible to vary the $\lambda^{a}$ and the AMOC reduction and sea-ice
response functions to tune the E-CCM even better to the CESM results. The
important result here is that the E-CCM reasonably captures the AMOC responses
from the CESM with only a few adjustment to the original model.

###  3.4 Transition Probabilities in the E-CCM

We define an on-to-off transition as a transition from the steady AMOC on
state to an AMOC state which strength is below a certain threshold (here taken
as 0 Sv for the no sea-ice case and 1 Sv for the sea-ice case). In addition,
an off-to-on transition is a transition from the steady AMOC off state to an
AMOC which strength is larger or equal to the strength of the steady AMOC on
state. We repeat the TAMS procedure (Lestang et al., 2018) 15 times and then
determine the mean transition probability $p$ (and its standard deviation) of
specific transitions within $T=100$ years, similar to Castellana et al.
(2019). Given this relatively short time interval (w.r.t. the dynamics of the
E-CCM) we only find typical F-type (fast) transitions (Castellana et al.,
2019). Full collapses (i.e., AMOC on $\rightarrow$ AMOC collapsed) or full
recoveries (i.e., AMOC collapsed $\rightarrow$ AMOC on), i.e. S-type
transitions (Castellana et al., 2019), take much longer time intervals (a few
centuries) and are not considered.

The CCM forced by freshwater noise ($f_{\sigma}=0.085$) for $E_{A}=0.25$ Sv
undergoes an on-to-off transition and then quickly returns to its AMOC on
state (Figures 6a). These AMOC transitions are reflected in the sign change of
the density difference between box n and box ts (Figures 6b). The AMOC state
is controlled by the sign in density difference and its magnitude determines
the transition probability. The magnitude of this density difference is larger
in the steady AMOC off state ($-$1.6 kg m-3) than the steady AMOC on state
(+0.7 kg m-3) and the larger the difference in absolute magnitude the lower
the transition probabilities (Figures 7a,b). In other words, realisations with
large density differences require relatively large freshwater noise
perturbation to make a transition. The transition probability under this
freshwater noise forcing is $p<10^{-9}$ for off-to-on transitions and is
substantially lower than for on-to-off transitions ($p=0.0025\pm 0.0001$).
This much lower transition probability also explains why there are no off-to-
on transitions for this particular noise realisation as more than $10^{9}$
noise realisations are required to expect at least 1 transition, whereas we
used only 5,000 realisations in Figures 6a,b. Lowering the freshwater forcing
towards the second bifurcation point ($E_{A}^{2}$ $\approx$0.06 Sv) would
favour off-to-on transitions while on the other hand the probability of on-to-
off transitions decreases (Figures 7c).

Figure 6: Transitions in the (E-)CCM for $f_{\sigma}=0.085$ and $E_{A}=0.25$
Sv for the AMOC strength (left column) and density difference between box n
and box ts (right column). The dashed lines indicate the statistical steady
AMOC on (red) and the AMOC off (blue) statistical steady state, where the
shading indicates the 95%-confidence level derived from 5,000 different noise
realisations (for each statistical steady state). A few random realisations
are displayed which (temporarily) undergo an on-to-off transition (red curves)
or an off-to-on transition (blue curves).

Under the same freshwater noise forcing for the E-CCM ($r_{q}=1$ (no sea-ice
effects), $\lambda^{a}=3.5\times 10^{-6}$ m s-1), we mainly find off-to-on
transitions (Figures 6c,d). In this set-up the system is much closer to its
second bifurcation point $E_{A}^{2}$; this is also reflected in the smaller
(absolute) density difference in the AMOC off state ($-$0.5 kg m-3) than in
the AMOC on state (+0.8 kg m-3). Under these freshwater noise forcing, the
transition probabilities are $p=7.9\times 10^{-7}\pm 0.9\times 10^{-7}$ (on-
to-off transitions) and $p=0.064\pm 0.002$ (off-to-on transitions). The
overall transition probabilities in the E-CCM with no sea-ice effects (Figures
7d,e,f) are higher than in the CCM (Figures 7a,b,c) when using the distance to
the bifurcation points as a reference for both model configurations. This is
because the magnitude of the freshwater noise forcing is given by
$f_{\sigma}E_{A}$ and, as the bifurcation points are shifted to higher values
of $E_{A}$ (Figure 3d) when temperature is included, it results in a larger
freshwater noise forcing allowing for higher transition probabilities.

The density differences remain unaltered ($<0.001\%$ difference w.r.t.
$r_{q}=1$) when including sea ice ($r_{q}^{\mathrm{min}}=0.02$) in the E-CCM.
Yet, there are now on-to-off transitions (Figures 6e,f). The density
difference of the two (random) realisations remains positive after an AMOC
collapse and the AMOC reduction factor (in the presence of sufficient sea ice)
is strong enough to limit their AMOC strength to almost 0 Sv. On relatively
short time scales ($<$ 100 years) the system undergoes an F-type transition
and remains in a weak AMOC state (i.e., AMOC strength slightly larger than 0
Sv). Note that this state is unstable for the imposed freshwater forcing
conditions and when the system further equilibrates, it makes the full
transition to the stable AMOC collapsed state. These transitions are somewhat
different compared to the other two model set-ups and we define on-to-off
transitions when $q_{N}^{\mathrm{ice}}<1$ Sv, and we label them as on-to-weak
transitions for the E-CCM. Sea ice suppresses the AMOC strength and this
strongly limits off-to-on transitions ($p<10^{-9}$) compared to on-to-weak
transitions ($p=0.0031\pm 0.0003$) under $f_{\sigma}=0.085$ and $E_{A}=0.25$
Sv. The much lower transition probabilities (Figures 7g,h,i) compared to the
other model set-ups clearly show that the AMOC weak/off state is much more
resilient than the AMOC on state. On the other hand, the on-to-weak
transitions are more likely when including sea ice (compare red curves in
Figures 7f and 7i for $E_{A}^{2}<E_{A}<E_{A}^{1}$) and this indicates that
sea-ice processes destabilise the AMOC on state.

Figure 7: Mean transition probabilities (within 100 years) in the (E-)CCM for
varying $f_{\sigma}$ and $E_{A}$, with $\Delta f_{\sigma}$ = 0.005 and $\Delta
E_{A}$ = 0.005 Sv. The on-to-off (or on-to-weak) transitions are shown in the
left column and the off(/weak)-to-on transitions in the middle column. For the
regions shaded in gray only one AMOC steady state exists. The dotted curve
indicates the freshwater noise forcing of 0.025 Sv (= $f_{\sigma}E_{A}$) and
the mean transition probabilities along this freshwater noise forcing are
shown in the right column, the shading indicates the minimum and maximum
transition probabilities derived from 15 independent TAMS realisations. The
inset in panel i shows the mean transition probabilities for the E-CCM for
$f_{\sigma}=0.3$ under varying $E_{A}$, and only for off/weak-to-on
transitions.

For $r_{q}^{\mathrm{min}}=0.02$ one stable AMOC weak state exists and this
state has both positive and negative $F_{\mathrm{ovS}}$ values (Figure 4).
Salinity perturbations are amplified under negative $F_{\mathrm{ovS}}$ (i.e.,
the salt-advection feedback) and this part of the AMOC weak state is
relatively unstable compared to the part with positive $F_{\mathrm{ovS}}$. The
quantity $\phi=\left|\frac{\partial\log p}{\partial E_{A}}\right|$ is indeed
larger for $F_{\mathrm{ovS}}<0$ (inset in Figure 7i) and for varying
$f_{\sigma}$, indicating that the salt-advection feedback destabilises this
part of the weak state where $F_{\mathrm{ovS}}<0$ and resulting in higher
transition probabilities.

Once the E-CCM undergoes an on-to-weak transition, it is also possible to
return back to the original AMOC on state. Figure 8 shows two realisations
which undergo such behaviour for $f_{\sigma}=0.25$ and $E_{A}=0.10$ Sv. Note
that we simulate the trajectories for 1,000 years in Figure 8 because the
model trajectory first needs to equilibrate to the AMOC weak state before the
second transition is possible. The strong transient responses in the
pycnocline and temperature and salinity for box ts dominate the dynamics of
the system as the AMOC strength is (partly) dependent on these oceanic
properties. The strong decline in the temperature and salinity (around model
year 200) is related to the deepening of the pycnocline, resulting in a volume
expansion of box ts and hence lowering of these quantities. The temperature
and salinity in box n (not shown) also decrease under an AMOC collapse and
remain fairly constant from model year 200 until AMOC recovery. Once the model
trajectory is sufficiently adjusted to the stable AMOC weak state (about 400
years after collapse), it may undergo again random transitions to the AMOC on
state. The transition probabilities from Figure 7 cannot (directly) be used to
determine the transition probability of undergoing two consecutive transitions
due to the setup of TAMS using fixed score functions (Castellana et al.,
2019).

Figure 8: Two random realisations undergoing transitions in the E-CCM
($\lambda^{a}=3.5\times 10^{-6}$ m s-1 and $r_{q}^{\mathrm{min}}=0.02$) for
$f_{\sigma}=0.25$ and $E_{A}=0.10$ Sv for various quantities. The dashed lines
indicate the steady AMOC on (red) and AMOC weak (black) states for each
presented quantity.

##  4 Summary and Discussion

In this study we analysed AMOC hysteresis behaviour in the CESM (van Westen
and Dijkstra, 2023) and branched off four additional simulations under
constant freshwater forcing to determine statistical steady states. We found
three different statistical steady states in the CESM: the AMOC on state, the
AMOC off state and the AMOC weak state. Sea-ice processes strongly influence
the AMOC strength and modify the AMOC hysteresis (van Westen and Dijkstra,
2023), giving rise to the AMOC weak state. This stable AMOC weak state is not
found in ocean-only models in which sea-ice processes are not resolved
(Dijkstra, 2007; Huisman et al., 2010). The existence of a broad multi-stable
regime under varying freshwater forcing suggests that noise-induced AMOC
tipping is possible in the CESM. However, we can not yet determine these
transition probabilities for the CESM due to computational limitations.

To determine the effects of the sea-ice processes on the AMOC transition
probabilities, we developed a novel stochastic ocean-sea-ice box model (the
E-CCM). The original 5-box model (Cimatoribus et al., 2014; Castellana et al.,
2019) was extended by dynamically solving for temperatures and parameterising
the effects of sea ice on the AMOC. When the effects of sea ice are strong
enough to suppress the AMOC strength, it modifies the AMOC hysteresis and
gives rise to the AMOC weak state. The deterministic box model can
qualitatively reproduce the CESM hysteresis behaviour, which is again an
example of the usefulness of conceptual models to understand AMOC dynamics
(Dijkstra, 2023). The existence of this weak state also illustrates the
crucial importance of sea-ice processes on the AMOC strength (van Westen and
Dijkstra, 2023). This is in agreement with the results in Lin et al. (2023),
who show in CMIP6 models that sea-ice effects play a major role (through their
effect on the Labrador Sea stratification) in the decline of the AMOC under
increased greenhouse gas concentrations.

The sea-ice parametrisation in the E-CCM only affects the AMOC strength and
does not consider brine rejection nor thermodynamics. Despite these
shortcomings, Hopf bifurcations exist in the E-CCM and this is consistent with
earlier works (Yang and Neelin, 1993; Brocker et al., 1990) where realistic
sea-ice AMOC interactions have been studied to understand natural AMOC
variability. Sea ice also plays an important role in Dansgaard-Oeschger events
(Peltier and Vettoretti, 2014; Vettoretti and Peltier, 2016; Boers et al.,
2018) and the recent AMOC hysteresis experiment conducted with the CESM
confirms this (van Westen and Dijkstra, 2023). Biases in the North Atlantic
sea-ice distribution influence the background AMOC strength and the transient
AMOC responses due to climate change (Lin et al., 2023) and can lead to rate-
induced tipping (Lohmann and Ditlevsen, 2021). These studies and the results
presented here suggest a prominent role for sea-ice processes in AMOC dynamics
and variability.

The stochastic E-CCM enables to compute the transition probabilities between
the different AMOC states under freshwater noise forcing using rare-event
algorithms (i.e., TAMS, Lestang et al. (2018)). The presence of sea ice
increases the probability of on-to-off and on-to-weak transitions and sea ice
strongly reduces the probability of weak-to-on and off-to-on transitions of
the AMOC. This indicates that the AMOC off and weak states, which are both
influenced by sea ice, are much more stable than the AMOC on state. Although
the quantitative results will be sensitive to the parameterisation of sea ice
and its effects on the AMOC in the E-CCM, they provide meaningful insights
into the effects of sea-ice processes in the AMOC hysteresis behaviour and
AMOC transition probabilities.

An interesting follow-up study is to analyse the AMOC transition behaviour,
and in particular the role of sea ice, under climate change. The suppressing
effects of sea ice on the AMOC become smaller as sea ice is less likely to
form under higher temperatures. As the AMOC weak state is strongly controlled
by sea ice, this state would (partly) disappear. This stabilises the AMOC as
transition probabilities from the AMOC on to AMOC weak state cease to exist
over an interval of freshwater forcing. On the other hand, climate change also
induces different temperature gradients which may destabilise the AMOC. A
combination of E-CCM analyses together with targeted CESM simulations can be
very useful for such a follow-up study.

###### Acknowledgements.

We thank Michael Kliphuis (IMAU, UU) for performing the additional CESM
simulations. The model simulations and the analysis of all the model output
was conducted on the Dutch National Supercomputer (Snellius) within NWO-SURF
project 17239. R.M.v.W., A.A.B. and H.A.D. are funded by the European Research
Council through the ERC-AdG project TAOC (PI: Dijkstra, project 101055096).
V.J.D. is funded by the European Union’s Horizon 2020 research and innovation
programme CriticalEarth under the Marie Sklodowska-Curie grant agreement
(project 956170). The authors declare no conflict of interest. The Python
codes for the E-CCM can be accessed through
https://doi.org/10.5281/zenodo.10554659, together with the CESM model output.
Supporting Information on the Idealised Model Simulations The CCM consists of
five equations for salinity (four dynamical and one to conserve salinity) and
one equation for the pycnocline. The equations are repeated here but we refer
to Castellana et al. (2019) where all the different terms are explained
($\theta$ is the Heaviside function) and input parameters are also provided:

$\displaystyle\frac{\mathrm{d}(V_{\mathrm{t}}S_{\mathrm{t}})}{\mathrm{d}t}$
$\displaystyle=q_{S}\left(\theta(q_{S})S_{\mathrm{ts}}+\theta(-q_{S})S_{\mathrm{t}}\right)+q_{U}S_{\mathrm{d}}-\theta(q_{N})q_{N}S_{\mathrm{t}}+r_{S}(S_{\mathrm{ts}}-S_{\mathrm{t}})+r_{N}(S_{\mathrm{n}}-S_{\mathrm{t}})+2E_{S}S_{0}$
(6a)
$\displaystyle\frac{\mathrm{d}(V_{\mathrm{ts}}S_{\mathrm{ts}})}{\mathrm{d}t}$
$\displaystyle=q_{Ek}S_{\mathrm{s}}-q_{e}S_{\mathrm{ts}}-q_{S}\left(\theta(q_{S})S_{\mathrm{ts}}+\theta(-q_{S})S_{\mathrm{t}}\right)+r_{S}(S_{\mathrm{t}}-S_{\mathrm{ts}})$
(6b)
$\displaystyle\frac{\mathrm{d}(V_{\mathrm{n}}S_{\mathrm{n}})}{\mathrm{d}t}$
$\displaystyle=\theta(q_{N})q_{N}\left(S_{\mathrm{t}}-S_{\mathrm{n}}\right)+r_{N}(S_{\mathrm{t}}-S_{\mathrm{n}})-(E_{S}+E_{A})S_{0}$
(6c)
$\displaystyle\frac{\mathrm{d}(V_{\mathrm{s}}S_{\mathrm{s}})}{\mathrm{d}t}$
$\displaystyle=q_{S}\left(\theta(q_{S})S_{\mathrm{d}}+\theta(-q_{S})S_{\mathrm{s}}\right)+q_{e}S_{\mathrm{ts}}-q_{Ek}S_{\mathrm{s}}-(E_{S}-E_{A})S_{0}$
(6d) $\displaystyle V_{\mathrm{d}}S_{\mathrm{d}}$
$\displaystyle=S_{0}V_{0}-V_{\mathrm{t}}S_{\mathrm{t}}-V_{\mathrm{ts}}S_{\mathrm{ts}}-V_{\mathrm{n}}S_{\mathrm{n}}-V_{\mathrm{s}}S_{\mathrm{s}}$
(6e)
$\displaystyle\left(A+\frac{L_{xA}L_{y}}{2}\right)\frac{\mathrm{d}D}{\mathrm{d}t}$
$\displaystyle=q_{U}+q_{Ek}-q_{e}-\theta(q_{N})q_{N}$ (6f)

The ocean heat content evolution equations (under the assumption of constant
heat capacity and reference density) are fairly similar to the salinity ones
and are given by:

$\displaystyle\frac{\mathrm{d}(V_{\mathrm{t}}T_{\mathrm{t}})}{\mathrm{d}t}$
$\displaystyle=q_{S}\left(\theta(q_{S})T_{\mathrm{ts}}+\theta(-q_{S})T_{\mathrm{t}}\right)+q_{U}T_{\mathrm{d}}-\theta(q_{N})q_{N}T_{\mathrm{t}}+r_{S}(T_{\mathrm{ts}}-T_{\mathrm{t}})+r_{N}(T_{\mathrm{n}}-T_{\mathrm{t}})-\lambda^{a}A_{\mathrm{t}}(T_{\mathrm{t}}-T_{\mathrm{t}}^{a})$
(7a)
$\displaystyle\frac{\mathrm{d}(V_{\mathrm{ts}}T_{\mathrm{ts}})}{\mathrm{d}t}$
$\displaystyle=q_{Ek}T_{\mathrm{s}}-q_{e}T_{\mathrm{ts}}-q_{S}\left(\theta(q_{S})T_{\mathrm{ts}}+\theta(-q_{S})T_{\mathrm{t}}\right)+r_{S}(T_{\mathrm{t}}-T_{\mathrm{ts}})-\lambda^{a}A_{\mathrm{ts}}(T_{\mathrm{ts}}-T_{\mathrm{ts}}^{a})$
(7b)
$\displaystyle\frac{\mathrm{d}(V_{\mathrm{n}}T_{\mathrm{n}})}{\mathrm{d}t}$
$\displaystyle=\theta(q_{N})q_{N}\left(T_{\mathrm{t}}-T_{\mathrm{n}}\right)+r_{N}(T_{\mathrm{t}}-T_{\mathrm{n}})-\lambda^{a}A_{\mathrm{n}}(T_{\mathrm{n}}-T_{\mathrm{n}}^{a})$
(7c)
$\displaystyle\frac{\mathrm{d}(V_{\mathrm{s}}T_{\mathrm{s}})}{\mathrm{d}t}$
$\displaystyle=q_{S}\left(\theta(q_{S})T_{\mathrm{d}}+\theta(-q_{S})T_{\mathrm{s}}\right)+q_{e}T_{\mathrm{ts}}-q_{Ek}T_{\mathrm{s}}-\lambda^{a}A_{\mathrm{s}}(T_{\mathrm{s}}-T_{\mathrm{s}}^{a})$
(7d)
$\displaystyle\frac{\mathrm{d}(V_{\mathrm{d}}T_{\mathrm{d}})}{\mathrm{d}t}$
$\displaystyle=\theta(q_{N})q_{N}T_{\mathrm{n}}-q_{S}\left(\theta(q_{S})T_{\mathrm{d}}+\theta(-q_{S})T_{\mathrm{s}}\right)-q_{U}T_{\mathrm{d}}$
(7e)

Using the standard input parameters from Castellana et al. (2019), we set the
ocean surface areas to $3\times 10^{13}$ m2 ($A_{\mathrm{s}}$), $1\times
10^{13}$ m2 ($A_{\mathrm{ts}}$), $10\times 10^{13}$ m2 ($A_{\mathrm{t}}$) and
$1\times 10^{13}$ m2 ($A_{\mathrm{n}}$). For more details how to integrate the
model we refer to the publicly-available Python software.

To initialise the E-CCM, we use the steady state solution from the CCM for an
arbitrary freshwater forcing $E_{A}$ (here 0.0644 Sv) in the AMOC on state.
The temperatures for box n and box ts are initialised at 5∘C and 10∘C
(Castellana et al., 2019), respectively, and the temperatures for the other
boxes can be arbitrarily chosen (here 0∘C). To remain at this steady state for
the E-CCM, we solve for each time step the equations
$\frac{\mathrm{d}(V_{\mathrm{ts}}T_{\mathrm{ts}})}{\mathrm{d}t}=0$ and
$\frac{\mathrm{d}(V_{\mathrm{n}}T_{\mathrm{n}})}{\mathrm{d}t}=0$ for
$T_{\mathrm{ts}}^{a}$ and $T_{\mathrm{n}}^{a}$, respectively. This guarantees
that the temperatures for box ts and box n remain the same over the entire
simulation and hence we remain at the same steady state as the one from
Castellana et al. (2019). To reduce the degrees of freedom, we assume the same
atmospheric temperatures for the tropical boxes
($T_{\mathrm{t}}^{a}=T_{\mathrm{ts}}^{a}$) and polar boxes
($T_{\mathrm{s}}^{a}=T_{\mathrm{n}}^{a}$). All the other temperatures are free
to evolve and we quickly reach a steady state for a given $\lambda^{a}$
(Figure 9a) and running the simulation further (2,000 years) does not change
the steady states (Figure 9b). The steady states of the deep ocean temperature
(not shown) are very close to $T_{\mathrm{n}}$ ($T_{\mathrm{s}}$) in the AMOC
on state (off state).

Box n receives relatively warm water from box t (in the AMOC on state), which
has to be balanced by a colder overhead atmosphere so that box n remains at a
constant temperature of 5∘C. For box ts, it is the exact opposite, this box
receives relatively cold water from box s (in the AMOC on state) and a warmer
overhead atmosphere is needed to balance this (negative) heat transport from
box s. The temperature of box s and box t are also influenced by the changing
atmospheric temperatures (they are coupled:
$T_{\mathrm{t}}^{a}=T_{\mathrm{ts}}^{a}$ and
$T_{\mathrm{s}}^{a}=T_{\mathrm{n}}^{a}$) and remain close to their overhead
atmospheric temperatures. Decreasing $\lambda^{a}$ to values below $3\times
10^{-6}$ m s-1 results in very cold (and unrealistic) atmospheric and oceanic
temperatures (Figure 9b). In the limit of $\lambda^{a}\rightarrow\infty$, the
atmospheric temperatures are 5∘C ($T_{\mathrm{n}}^{a}$) and 10∘C
($T_{\mathrm{ts}}^{a}$) and temperature anomalies in box n and box ts decay
immediately, resulting in constant temperatures of 5∘C and 10∘C, respectively.
This gives the same results as in the CCM. After determining the steady states
for a given $\lambda^{a}$ (Figure 9b), we fix the atmospheric temperatures and
use those as input parameters for the E-CCM. All ocean temperatures are now
free to evolve under a changing freshwater forcing $E_{A}$ and the steady
state is identical to the one from Castellana et al. (2019) for the reference
$E_{A}$ value (here 0.0644 Sv).

Figure 9: (a): Evolution of the ocean (solid) and atmosphere (dashed)
temperatures over the first 100 years for $\lambda^{a}=5\times 10^{-6}$ m s-1,
where $T_{\mathrm{n}}=5^{\circ}$C and $T_{\mathrm{ts}}=10^{\circ}$C are fixed.
(b): The steady states (after 2,000 years of integration) temperatures for
varying $\lambda^{a}$. Figure 10: (a): Sea-ice response function as a function
of the temperature of box n. (b): The AMOC reduction factor as a function of
the sea-ice fraction of box n for a given $r_{q}^{\mathrm{min}}$.

## References

  * Armstrong McKay et al. (2022) Armstrong McKay, D. I., and Coauthors, 2022: Exceeding 1.5 C global warming could trigger multiple climate tipping points. Science, 377 (6611), eabn7950.
  * Ashwin et al. (2012) Ashwin, P., S. Wieczorek, R. Vitolo, and P. Cox, 2012: Tipping points in open systems: bifurcation, noise-induced and rate-dependent examples in the climate system. Philosophical Transactions of the Royal Society A, 370, 1166–1184.
  * Baars et al. (2021) Baars, S., D. Castellana, F. Wubs, and H. Dijkstra, 2021: Application of adaptive multilevel splitting to high-dimensional dynamical systems. Journal of Computational Physics, 424, 109 876, 10.1016/j.jcp.2020.109876, 2011.05745.
  * Baatsen et al. (2020) Baatsen, M., A. S. Von Der Heydt, M. Huber, M. A. Kliphuis, P. K. Bijl, A. Sluijs, and H. A. Dijkstra, 2020: The middle to late Eocene greenhouse climate modelled using the CESM 1.0.5. Climate of the Past, 16 (6), 2573–2597.
  * Boers (2021) Boers, N., 2021: Observation-based early-warning signals for a collapse of the Atlantic Meridional Overturning Circulation. Nature Climate Change, 11 (8), 680–688, 10.1038/s41558-021-01097-4.
  * Boers et al. (2018) Boers, N., M. Ghil, D. R. P. o. the, and 2018, 2018: Ocean circulation, ice shelf, and sea ice interactions explain Dansgaard–Oeschger cycles. Proceedings National Acad Sciences, 10.1073/pnas.1802573115/-/dcsupplemental.
  * Brocker et al. (1990) Brocker, W., G. Bond, M. Klas, G. Bonani, and W. Wolfli, 1990: A salt oscillator in the glacial atlantic? 1. The concept: Paleoceanography, 5, 469–477.
  * Caesar et al. (2021) Caesar, L., G. D. McCarthy, D. Thornalley, N. Cahill, and S. Rahmstorf, 2021: Current Atlantic meridional overturning circulation weakest in last millennium. Nature Geoscience, 14 (3), 118–120.
  * Caesar et al. (2018) Caesar, L., S. Rahmstorf, A. Robinson, G. Feulner, and V. Saba, 2018: Observed fingerprint of a weakening Atlantic Ocean overturning circulation. Nature, 556 (7700), 191–196.
  * Castellana et al. (2019) Castellana, D., S. Baars, F. W. Wubs, and H. A. Dijkstra, 2019: Transition probabilities of noise-induced transitions of the Atlantic Ocean circulation. Scientific Reports, 9 (1), 20 284.
  * Cessi (1994) Cessi, P., 1994: A simple box model of stochastically forced thermohaline flow. Journal of Physical Oceanography, 24, 1911–1920.
  * Cimatoribus et al. (2014) Cimatoribus, A. A., S. S. Drijfhout, and H. A. Dijkstra, 2014: Meridional overturning circulation: Stability and ocean feedbacks in a box model. Climate dynamics, 42, 311–328.
  * Dijkstra (2007) Dijkstra, H. A., 2007: Characterization of the multiple equilibria regime in a global ocean model. Tellus A: Dynamic Meteorology and Oceanography, 59 (5), 695–705.
  * Dijkstra (2023) Dijkstra, H. A., 2023: The role of conceptual models in climate research. Physica D: Nonlinear Phenomena, 133984.
  * Ditlevsen and Ditlevsen (2023) Ditlevsen, P., and S. Ditlevsen, 2023: Warning of a forthcoming collapse of the Atlantic meridional overturning circulation. Nature Communications, 14 (1), 4254.
  * Doedel et al. (2007) Doedel, E. J., R. C. Paffenroth, A. C. Champneys, T. F. Fairgrieve, Y. A. Kuznetsov, B. E. Oldeman, B. Sandstede, and X. J. Wang, 2007: AUTO-07p: Continuation and Bifurcation Software for Ordinary Differential Equations.
  * Doedel et al. (2021) Doedel, E. J., R. C. Paffenroth, A. C. Champneys, T. F. Fairgrieve, Y. A. Kuznetsov, B. E. Oldeman, B. Sandstede, and X. J. Wang, 2021: auto-07p. URL https://github.com/auto-07p/auto-07p.
  * Garzoli et al. (2013) Garzoli, S. L., M. O. Baringer, S. Dong, R. C. Perez, and Q. Yao, 2013: South Atlantic meridional fluxes. Deep Sea Research Part I: Oceanographic Research Papers, 71, 21–32.
  * Hawkins et al. (2011) Hawkins, E., R. S. Smith, L. C. Allison, J. M. Gregory, T. J. Woollings, H. Pohlmann, and B. De Cuevas, 2011: Bistability of the Atlantic overturning circulation in a global climate model and links to ocean freshwater transport. Geophysical Research Letters, 38 (10), L10 605.
  * Hu et al. (2012) Hu, A., and Coauthors, 2012: Role of the Bering Strait on the hysteresis of the ocean conveyor belt circulation and glacial climate stability. Proceedings of the National Academy of Sciences, 109 (17), 6417–6422, 10.1073/pnas.1116014109.
  * Huisman et al. (2010) Huisman, S. E., M. Den Toom, H. A. Dijkstra, and S. Drijfhout, 2010: An indicator of the multiple equilibria regime of the atlantic meridional overturning circulation. Journal of Physical Oceanography, 40 (3), 551–567.
  * Jackson et al. (2022) Jackson, L. C., A. Biastoch, M. W. Buckley, D. G. Desbruyères, E. Frajka-Williams, B. Moat, and J. Robson, 2022: The evolution of the North Atlantic meridional overturning circulation since 1980. Nature Reviews Earth & Environment, 3 (4), 241–254.
  * Jacques-Dumas et al. (2023) Jacques-Dumas, V., R. M. van Westen, F. Bouchet, and H. A. Dijkstra, 2023: Data-driven methods to estimate the committor function in conceptual ocean models. Nonlinear Processes in Geophysics, 30 (2), 195–216.
  * Lestang et al. (2018) Lestang, T., F. Ragone, C.-E. Bréhier, C. Herbert, and F. Bouchet, 2018: Computing return times or return periods with rare event algorithms. Journal of Statistical Mechanics: Theory and Experiment, 2018 (4), 043 213.
  * Lin et al. (2023) Lin, Y.-J., B. E. Rose, and Y.-T. Hwang, 2023: Mean state AMOC affects AMOC weakening through subsurface warming in the Labrador Sea. Journal of Climate, 36 (12), 3895–3915.
  * Lohmann and Ditlevsen (2021) Lohmann, J., and P. D. Ditlevsen, 2021: Risk of tipping the overturning circulation due to increasing rates of ice melt. Proceedings of the National Academy of Sciences, 118 (9), e2017989 118, 10.1073/pnas.2017989118.
  * Marotzke (2000) Marotzke, J., 2000: Abrupt climate change and thermohaline circulation: Mechanisms and Predictability. Proc. Natl. Acad. Sci., 97, 1347–1350.
  * Mecking et al. (2017) Mecking, J., S. Drijfhout, L. Jackson, and M. Andrews, 2017: The effect of model bias on Atlantic freshwater transport and implications for AMOC bi-stability. Tellus A: Dynamic Meteorology and Oceanography, 69 (1), 1299 910.
  * Mehling et al. (2023) Mehling, O., R. Börner, and V. Lucarini, 2023: Limits to predictability of the asymptotic state of the Atlantic Meridional Overturning Circulation in a conceptual climate model. arXiv preprint arXiv:2308.16251.
  * Orihuela-Pinto et al. (2022) Orihuela-Pinto, B., M. H. England, and A. S. Taschetto, 2022: Interbasin and interhemispheric impacts of a collapsed atlantic overturning circulation. Nature Climate Change, 12 (6), 558–565.
  * Peltier and Vettoretti (2014) Peltier, W. R., and G. Vettoretti, 2014: Dansgaard-Oeschger oscillations predicted in a comprehensive model of glacial climate: A “kicked” salt oscillator in the Atlantic. Geophysical Research Letters, 41 (20), 7306–7313.
  * Rahmstorf et al. (2005) Rahmstorf, S., and Coauthors, 2005: Thermohaline circulation hysteresis: a model intercomparison. Geophysical Research Letters, L23605, 1–5.
  * Smeed et al. (2018) Smeed, D. A., and Coauthors, 2018: The North Atlantic Ocean is in a state of reduced overturning. Geophysical Research Letters, 45 (3), 1527–1533.
  * Stommel (1961) Stommel, H., 1961: Thermohaline convection with two stable regimes of flow. Tellus, 13 (2), 224–230.
  * van Westen and Dijkstra (2023) van Westen, R. M., and H. A. Dijkstra, 2023: Asymmetry of AMOC Hysteresis in a State-of-the-Art Global Climate Model. Geophysical Research Letters, 50 (22), e2023GL106 088.
  * van Westen et al. (2024) van Westen, R. M., M. Kliphuis, and H. A. Dijkstra, 2024: Physics-Based Early Warning Signal shows AMOC is on Tipping Course. Science Advances, in press,, 1–29.
  * Vettoretti and Peltier (2016) Vettoretti, G., and W. R. Peltier, 2016: Thermohaline instability and the formation of glacial North Atlantic super polynyas at the onset of Dansgaard-Oeschger warming events. Geophysical Research Letters, 43 (10), 5336–5344.
  * Weijer et al. (2020) Weijer, W., W. Cheng, O. A. Garuba, A. Hu, and B. T. Nadiga, 2020: CMIP6 Models Predict Significant 21st Century Decline of the Atlantic Meridional Overturning Circulation. Geophysical Research Letters, 47 (12), e2019GL08 607, 10.1029/2019gl086075.
  * Worthington et al. (2021) Worthington, E. L., B. I. Moat, D. A. Smeed, J. V. Mecking, R. Marsh, and G. D. McCarthy, 2021: A 30-year reconstruction of the Atlantic meridional overturning circulation shows no decline. Ocean Science, 17 (1), 285–299.
  * Yang and Neelin (1993) Yang, J., and J. D. Neelin, 1993: Sea-ice interaction with the thermohaline circulation. Geophys.Res. Letters, 20, 217 – 220.