Characterisation of Atlantic meridional overturning hysteresis using Langevin dynamics

Hysteresis diagrams of the Atlantic meridional overturning circulation (AMOC) under freshwater forcing from climate models of intermediate complexity are fitted to a simple model based on the Langevin equation. A total of six parameters are sufficient to quantitatively describe the collapses seen in these simulations. Reversing the freshwater forcing results in asymmetric behaviour that is less well captured and appears to require a more complicated model. The Langevin model allows for comparison between models that display an AMOC collapse. Differences between the climate models studied here are mainly due to the strength of the stable AMOC and the strength of the response to a freshwater forcing.

1 Introduction

The Atlantic meridional overturning circulation (AMOC) is an important circulation in the Atlantic Ocean. It is also an important part of the climate system overall due to the heat it transports from the South Atlantic to the North Atlantic (Ganachaud and Wunsch, 2000; Vellinga and Wood, 2002). The AMOC therefore has a substantial influence on the (western) European climate, and a weakening of the AMOC might cause changes in the European climate and weather. The AMOC has also been identified as one of Earth's “tipping elements”, where a rapid change on markedly faster timescales could take place in the (near) future (Lenton et al., 2008). The AMOC is partly buoyancy-driven by the deep water formations in the North Atlantic subpolar gyre, which produces the North Atlantic Deep Water (NADW) (e.g. Rahmstorf, 2000). The AMOC might be bistable in nature, which means it admits an “off” state, with little or no transport from north to south, as a counterpart to its current “on” state (Broecker et al., 1985).

Palaeoclimate records of the last glacial period show a rapid switching of temperature, which might be associated with the presence or absence of a vigorous AMOC as it exists today (Dansgaard et al., 1993). The possibility of a bistable AMOC being the cause of these rapid changes has been noted (Broecker et al., 1990). With the current climate warming rapidly, the stability of the AMOC is of particular interest (Collins et al., 2013), and climate modelling projections indicate that AMOC strength will decrease under an increase of CO₂. Recent measurements show the AMOC has decreased in strength (Smeed et al., 2018). An understanding of the possibly bistable nature of the AMOC is therefore relevant to understand the consequences of climate change. See Weijer et al. (2019) for a review on AMOC bistability.

The Langevin equation has been posited before as being suitable to capture the essential dynamics of an AMOC collapse (Ditlevsen and Johnsen, 2010; Berglund and Gentz, 2002). It has also been used elsewhere as the basis for describing the dynamics of climate subsystems (Kwasniok and Lohmann, 2009; Livina et al., 2010) and the AMOC in particular (Kleinen et al., 2003; Held and Kleinen, 2004). A fourth-order potential function is used in Ditlevsen and Johnsen (2010); Berglund and Gentz (2002) because it is the minimum required for having three distinct solutions (double wells). This potential function has two parameters, which are presumed to be functions of the freshwater forcing. Variation in the freshwater forcing is assumed to directly drive changes in AMOC strength by changing the potential function in the Langevin equation. Although the hysteresis loops of the AMOC include both a collapse and a resurgence point, we will only attempt to model the collapse from the stable “on” branch to the stable “off” branch.

Though the Langevin equation has played a role in the conceptual picture of bistability and tipping points in the climate, it has not been used to actually fit the parameters to a (simulated) AMOC collapse. Here, we attempt to construct a simple model based on the Langevin equation and fit its dynamics to salt-advection-driven collapse trajectories of the AMOC seen in climate models (Rahmstorf et al., 2005). The result is a set of parameters that quantitatively describe the AMOC collapse process. This derived model defines a low-dimensional manifold that captures the essential AMOC collapse characteristics. To the extent that the low-dimensional model is successful in capturing the more complex model, this method could also be used to predict the parameter range where in a model a collapse would occur. At present, however, it is intended to provide a characterisation of the collapse that will allow comparison between climate models.

Section 2 sketches the theoretical background of the Langevin equation and of the salt-advection mechanism. In Sect. 3 we fit the proposed Langevin model to the AMOC collapse trajectories seen in a set of climate models of intermediate complexity (EMICs) taken from Rahmstorf et al. (2005). We end with a discussion and conclusions in Sect. 4.

2 The Langevin model

An increase in surface air temperatures or an increased surface freshwater flux by changes in precipitation minus evaporation will decrease the buoyancy in the shallow layer of the deep water formation regions in the North Atlantic subpolar gyre. The deep water formation is reduced, and the southward meridional flow is also reduced. In principle, this mechanism can reduce the AMOC to zero gradually if fully buoyancy-driven. A salt-advection feedback mechanism that leads to a bimodal AMOC was proposed by Stommel (1961). In this mechanism, salinity anomalies in the North Atlantic are amplified by the overturning flow, which in turn controls the North Atlantic salinity. Positive anomalies are strengthened and negative anomalies weakened; this results in a positive feedback between the salinity anomalies and the overturning. Bistability, consisting of a strong and a weak AMOC state, and possible abrupt transitions result from this process.

Figure 1 Example bifurcation diagram of the AMOC (Ψ) in response to a control variable μ. The red branch is the on state (upper), and blue is the off state (lower). The upper branch deforms when closer to the bifurcation points that are connected though the repeller (dashed line). The two bifurcations points are indicated as μ₊ (collapse point) and μ₋ (resurgence point). The top ± symbols indicate a unimodal (+) or bimodal (−) regime.

Figure 1 shows a conceptual picture of the two stable AMOC (index) states. The AMOC is a scalar variable obtained by integrating the overturning transport and selecting its maximum value (typically located in the subtropical North Atlantic). In red, the upper branch is drawn up to the collapse point where a bifurcation occurs. The real AMOC in the current climate moves along this branch from the left to the right, towards its (assumed) collapse point. The branch in blue is the counterpart of the upper branch and represents the off state of the AMOC and ends in another bifurcation point to the left where the AMOC jumps back to full strength. The dashed line (repeller) separates the two basins of attraction associated with the two stable branches (attractors). At the bifurcation point one of the two basins of attraction vanishes and a qualitative change takes place in the potential function (the number of solutions for a given value of the freshwater forcing μ goes from 3 to 1).

Below we will derive a model based on the Langevin equation that captures the essential dynamics of a bimodal AMOC under a freshwater forcing μ.

2.1 Multiple stable AMOC states

The conceptual picture of the AMOC being a zero-dimensional variable that is driven by stochastic forces trapped in a potential is similar to that of a particle's motion described by Langevin dynamics (Lemons et al., 1908). The Langevin equation (Gardiner, 2004; Ditlevsen and Johnsen, 2010) is as follows:

$$\begin{array}{}\text{(1)}& \dot{x}=-{\partial }_x{U}_{\mathit{\mu }}\left(x\right)+\mathit{\sigma }\mathit{\eta }.\end{array}$$

It describes the position of a noise-driven particle (x) trapped in a potential function U. The stochastic term is a white noise process (η) scaled with an intensity parameter σ. At first we will ignore the stochastic nature of the AMOC collapse process and focus on the deterministic behaviour.

The double well potential seen in Fig. 1 has been extensively studied and applied in a quantitative way as well. However, to our knowledge it has not been quantitatively applied to AMOC hysteresis using the Langevin equation in complex numerical climate models before.

AMOC bistability has, however, been studied quantitatively in, e.g. Boulton et al. (2014), using transient runs. In Poston and Stewart (1978) an extensive treatment is given why, in addition to a scaling and shifting, only two parameters are sufficient to describe the bistability. More precisely, the third-order term and the fourth-order coefficient can be eliminated. The two remaining coefficients in the polynomial describe the critical behaviour, not just locally near the critical points but also for the entire trajectory under a suitable transformation. A direct consequence is that only partial information, in the form of a piece of the trajectory, should suffice to describe the entire trajectory (the full hysteresis loop).

The potential function takes the following form (Gardiner, 2004; Ditlevsen and Johnsen, 2010):

$$\begin{array}{}\text{(2)}& -U\left(x\right)=-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{x}^{\mathrm{4}}+{\displaystyle \frac{\mathit{\beta }}{\mathrm{2}}}{x}^{\mathrm{2}}+\mathit{\alpha }x.\end{array}$$

The two parameters α and β are functions of the freshwater forcing μ. The AMOC state variable Ψ requires an affine transformation (Cobb, 1980),

$$\begin{array}{ll}{\displaystyle }\mathit{\alpha }& {\displaystyle }=\mathit{\alpha }\left(\mathit{\mu }\right)\\ {\displaystyle }\mathit{\beta }& {\displaystyle }=\mathit{\beta }\left(\mathit{\mu }\right)\\ {\displaystyle }x& {\displaystyle }=(\mathrm{\Psi }-\mathit{\lambda })/\mathit{\nu }.\end{array}$$

To fit the model trajectories we need to find expressions for α and β and suitable values for the transformation parameters λ and ν. In the literature α is referred to as the normal factor, and β is referred to as the splitting factor (Poston and Stewart, 1978). In the bifurcation diagram the value of ν is approximately the distance in Ψ between the bifurcation point on the top branch to the bifurcation point on the lower branch. Similarly, the value of λ is approximately the Ψ value between the bifurcation points at μ_±. The transformation uses λ to shift the trajectory and ν to scale it. Below we describe the potential visually and state additional constraints that follow from the demand that the freshwater forcing is the only variable that determines the dynamical behaviour.

2.2 Potential description

In Fig. 2 an overview of the qualitatively different forms of potential are shown (−U(x), right column) together with their derivative functions ($-{\partial }_{x}U$, left column). Dots indicate the location of critical points and are related to the number of wells in the potential. The top panels show the typical bimodal form (I) with two stable states and one unstable state in the middle. Below these are the three possible unimodal states (E). These occur for forcing values to the left of μ₋ and to the right of μ₊. The panels B₁ and B₂ are the sub-manifolds that separates the unimodal regime from the bimodal regime. These two meet in the cusp point P, as shown in the bottom panels. See Poston and Stewart (1978) for further details.

Figure 2 Sample potentials (right) and their derivatives (left) for (top to bottom) the three possible varieties of bimodal state (I), three types of unimodal state (E), the two pathological cases where D=0 (B₁ and B₂), and the cusp catastrophe point (P). Dots indicate the critical points. Scaling is not uniform between panels. Note the choice of negative sign of the potential U.

In Fig. 3 the stability diagram is shown where the areas indicated are those with qualitatively different behaviour seen in Fig. 2. See Poston and Stewart (1978) for similar diagrams. The cusp point P is the singular point where no proper solution can exist because only the trivial solution (all parameters are valued 0) is allowed here (both bifurcation points μ_± and AMOC strength are at zero). The two parameters are α and β and are the two coefficients in the potential function. Their values change because of their dependency on the forcing value (μ).

Figure 3 Discriminant determining the stability and number of critical points. The splitting factor β and normal factor α describe the stability diagram. The bimodal regime (I) is separated from the unimodal regime (E) by two lines (B_1, 2) that meet at point P.

Our aim is to arrive at a description that matches a series of μ values across the stability diagram. The two parameters α and β are independent but can be parameterised by other variables that map them to observations. If parameterised by a single variable, the values of (α, β) across the stability surface are a one-dimensional subset, as suggested by the AMOC index. On one side of the cusp point, along the splitting axis (β), only a unimodal regime exists, while on the other side two regimes exist with the modes at relative distances apart.

2.3 Constraints

With a varying α there exist an interval between two critical points (α_±) in between which the distribution is bimodal and unimodal outside that interval. Because the AMOC trajectory is one-dimensional and μ is also one-dimensional, there must be a relation between α and β that reduces dimensionality from two dimensions to one dimension. When passing through the critical point α₊, the number of potential wells goes from two to one. Similarly, moving through α₋ changes the number of wells from one to two (for given μ_±). The two critical points of ∂_xU and μ_± can be found analytically for μ_± real and degenerate solutions. It can be shown (Birkhoff and Mac Lane, 1970, p. 106) that the discriminant $D=\mathrm{27}{\mathit{\alpha }}^{\mathrm{2}}-\mathrm{4}{\mathit{\beta }}^{\mathrm{3}}=\mathrm{0}$ (i.e. real solutions) needs to be solved for α to obtain the two critical solutions that relate α and β. It is at these solutions that the number of critical points changes at forcing values μ_±. When D<0, there are three distinct real solutions, which correspond to the bimodal regime; when D>0, there is only one distinct real solution, which corresponds to the unimodal regime. When any two of the roots are the same, the number of extrema goes from three to two (or one if all are the same) and the solutions become degenerate (this occurs at B_1, 2 in Fig. 3).

Solving for α gives two solutions that are the critical values as functions of β,

$${\mathit{\alpha }}_{\pm }=\pm {\displaystyle \frac{\mathrm{2}\sqrt{\mathrm{3}}}{\mathrm{9}}}{\left(\mathit{\beta }\right)}^{\mathrm{3}/\mathrm{2}}\text{or}{\mathit{\alpha }}_{\pm }=\mp {\displaystyle \frac{\mathrm{2}\sqrt{\mathrm{3}}}{\mathrm{9}}}{\left(\mathit{\beta }\right)}^{\mathrm{3}/\mathrm{2}},$$

with β≥0 for real solutions. The points α_± correspond to where the lines B_1, 2 in Fig. 3 are passed when moving across the stability surface.

For ${\mathit{\alpha }}_{+}<\mathrm{0}$ $-U\left(\mathrm{1}\right)<\mathrm{0}$, this corresponds with the AMOC undergoing a collapse at μ₊ from an on state to an off state, and the correct choice of sign is

$$\begin{array}{}\text{(3)}& {\mathit{\alpha }}_{\pm }=\mp {\displaystyle \frac{\mathrm{2}\sqrt{\mathrm{3}}}{\mathrm{9}}}{\left({\mathit{\beta }}_{\pm }\right)}^{\mathrm{3}/\mathrm{2}},\end{array}$$

with α_± and β_± being the values corresponding to μ_±. Changing μ in the bifurcation diagram corresponds to moving from curve B₂ to curve B₁, and Eq. (3) relates the two stability parameters α and β at the two critical forcing values μ_±.

2.3.1 Linear functions α,β

The value of β does not need to be fixed (to α_± and in general there is a corresponding β_± at the respective critical points. We assume linear functions for α and β,

$$\begin{array}{l}{\displaystyle }\mathit{\alpha }\left(\mathit{\mu }\right)={\mathit{\alpha }}_{\mathrm{0}}+\mathit{\mu }\mathit{\delta }\mathit{\alpha },\\ {\displaystyle }\mathit{\beta }\left(\mathit{\mu }\right)={\mathit{\beta }}_{\mathrm{0}}+\mathit{\mu }\mathit{\delta }\mathit{\beta },\end{array}$$

reducing the dependency to these four parameters. Linear functions are the simplest non-trivial dependencies, while adding non-linear parameters introduces further unknowns, making this the most parsimonious parameterisation that captures the first-order behaviour. In addition, intuitively we can understand the pair (δα, δβ) as the angle under which the system moves to the bifurcation point (B_1, 2) in Fig. 3), which locally only requires the values of α and β up to first order. From this parameterisation we can determine the offset α₀ and rate δα in terms of β₀ and δβ,

$$\begin{array}{rl}& {\mathit{\alpha }}_{+}={\mathit{\alpha }}_{\mathrm{0}}+{\mathit{\mu }}_{+}\mathit{\delta }\mathit{\alpha }=-{\displaystyle \frac{\mathrm{2}\sqrt{\mathrm{3}}}{\mathrm{9}}}{\left({\mathit{\beta }}_{+}\right)}^{\mathrm{3}/\mathrm{2}}\text{and}\\ & {\mathit{\alpha }}_{-}={\mathit{\alpha }}_{\mathrm{0}}+{\mathit{\mu }}_{-}\mathit{\delta }\mathit{\alpha }=+{\displaystyle \frac{\mathrm{2}\sqrt{\mathrm{3}}}{\mathrm{9}}}{\left({\mathit{\beta }}_{-}\right)}^{\mathrm{3}/\mathrm{2}},\end{array}$$

gives

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathit{\delta }\mathit{\alpha }=-{\displaystyle \frac{\mathrm{2}\sqrt{\mathrm{3}}}{\mathrm{9}}}{\displaystyle \frac{{\left({\mathit{\beta }}_{+}\right)}^{\mathrm{3}/\mathrm{2}}+{\left({\mathit{\beta }}_{-}\right)}^{\mathrm{3}/\mathrm{2}}}{{\mathit{\mu }}_{+}-{\mathit{\mu }}_{-}}},\text{(5)}& {\displaystyle }\begin{array}{rl}{\mathit{\alpha }}_{\mathrm{0}}=\mathit{\alpha }(\mathit{\mu }=\mathrm{0})& ={\displaystyle \frac{\sqrt{\mathrm{3}}}{\mathrm{9}}}\left[-{\left({\mathit{\beta }}_{+}\right)}^{\mathrm{3}/\mathrm{2}}+{\left({\mathit{\beta }}_{-}\right)}^{\mathrm{3}/\mathrm{2}}\right]\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\delta }\mathit{\alpha }\left({\mathit{\mu }}_{+}+{\mathit{\mu }}_{-}\right).\end{array}\end{array}$$

This constrains the values of α, leaving only β as a free variable, which is then parameterised by β₀ and δβ. Note that only solutions with ${\mathit{\beta }}_{\pm }>\mathrm{0}$ are valid. In addition, values for β₀ and δβ₀ that result in crossing B₂ in another point besides β₋ are unsuitable. The curves B_1, 2 are each intersected by a straight line in at most two points, and we require intersection at a single point only.

2.4 Stochastic interpretation

With the deterministic framework in place, the stochastic nature can be reintroduced. The potential function can be replaced by a distribution that is the stationary distribution in the asymptotic limit (i.e. the long-term behaviour of repeated sampling of the hysteresis loop). The potential (a fourth-order polynomial) gives the following probability distribution (Cobb, 1978):

$$\begin{array}{}\text{(6)}& P(x,\mathit{\alpha },\mathit{\beta })=C{e}^{-\mathrm{2}/{\mathit{\sigma }}^{\mathrm{2}}U\left(x\right)}=C{e}^{\mathrm{2}/{\mathit{\sigma }}^{\mathrm{2}}(-\mathrm{1}/\mathrm{4}{x}^{\mathrm{4}}+\mathit{\beta }/\mathrm{2}{x}^{\mathrm{2}}+\mathit{\alpha }x)}.\end{array}$$

The factor $C=C(\mathit{\alpha },\mathit{\beta })$ does not have a (known) analytical expression for the general case but can be computed numerically (and can therefore used as a factor in the likelihood function in the next section). This can be done accurately with an adaptive quadrature method (Piessens et al., 2012), though it suffers from numerical limitations. The value of σ is a measure of intrinsic variation in the AMOC. Note that σ is a measure of additive noise (because we assume that σ is not dependent on μ), and other choices, such as multiplicative noise, can be made (Das and Kantz, 2020). See Gardiner (2004) for a derivation of this distribution using the Fokker–Planck equation, from which the Langevin equation can also be derived. In addition, note that $\mathit{\sigma }\to \mathit{\sigma }/\mathit{\nu }$ because of the scaling with ν we introduced in Sect. 2.1.

Figure 4 Example trajectory with corresponding distribution. Parameterised by λ=15, ν=20, σ=0.12ν, ${\mathit{\mu }}_{+}=\mathrm{0.2}$, ${\mathit{\mu }}_{-}=\mathrm{0}$, β₀=0.2, and δβ=0; α₀ and δα follow the constraints in Eqs. (4) and (5). The distribution of one of the attractor branches (red: on state; blue: off state) deforms when closer to the bifurcation points that are connected though the repeller that forms the trench of the distribution (dashed line). Top ± symbols indicate a unimodal (+) or bimodal (−) regime based on the discriminant value (D). The value of σ is relatively large and is chosen for clarity. The purple lines indicate the (fixed) positions of the bifurcation points.

An example bifurcation diagram with corresponding distribution is shown in Fig. 4. The purple lines indicate the (fixed) positions of the bifurcation points. The dashed grey line marks the positions of the unstable solution (repeller) in between the two attractor branches that separates the two basins of attraction. Note that the bifurcation points are extremal in the sense that no bimodality can exist beyond them. With the trajectories being noisy and driven along the attractor, there is (always) some probability of a “noise-induced” transition. The state shifts from one basin of attraction to the other, crossing the repeller, and the AMOC rapidly moves from one attractor to the other. For this reason, the bimodality region might be larger than is apparent from a particular sample AMOC trajectory. A larger noise level (as seen in AMOC observations Smeed et al., 2018) would increase the likelihood of a collapse before the AMOC reaches the bifurcation point.

Figure 5 (a) Distributions from the exponential family (Eq. 6) where the parameter β is kept at a fixed value and α is varied. The distribution transforms from unimodal (back), to bimodal (middle), to a different unimodal distribution (front). The bimodal states have a larger and a smaller mode, depending on the position within the bimodal regime. The relative strength between modes depends on σ. (b) Distributions from the exponential family (Eq. 6) where the parameter α is kept at a fixed value and β is varied. A broad unimodal state (at the back) splits into distinct bimodal states (to the front). In the middle a critical point exists, called the cusp (point P in Fig. 3), where the split occurs.

The distributions in Fig. 5 show that qualitatively distinct behaviour occurs when α or β are varied. For both parameters, a change from a unimodal to a bimodal distribution can be seen. Each distinct shape of the distribution can be identified with one of the potential functions in Fig. 2. In principle, a change in only one of the two structural parameters (α and β) can move the distribution between unimodal and bimodal forms.

We are now in a position to apply the above to collapse trajectories from climate models.

3 AMOC collapse parameter estimation

We describe how to find an optimal solution under the framework described in the previous section. Using a Bayesian optimisation procedure, estimated values of β₀ and δβ can be found, together with the scaling parameters ν and λ. We will also estimate the values for μ_±, resulting in a six-parameter list that describes (the upper branch) of an AMOC collapse.

The parameters β₀ and δβ are independent of each other but need to cross the curves B_1, 2 in Fig. 3) to match the corresponding values for μ_±. This constraint is satisfied by the resulting values for α₀ and δα. This can still lead to solution candidates that are not suitable for the collapse trajectories and are eliminated in the sampling process below. The scaling parameters are not fully independent because λ<ν (the offset cannot exceed the scaling) and knowing where the upper and lower branches are located already gives a rough estimate.

3.1 Parameter estimation

Cobb (1978) was able to fit the distribution in Eq. (6) using optimisation techniques (which were numerically unstable and not very flexible). Though the estimates for the scaling parameters λ and ν can be quite good with this approach, estimating the trajectory parameters β₀ and δβ requires a more flexible method. Knowing which distribution to use, we can estimate the posterior probability distribution of the parameters given the data Ψ(μ),

$$P(\mathit{\nu },\mathit{\lambda },{\mathit{\beta }}_{\mathrm{0}},\mathit{\delta }\mathit{\beta },{\mathit{\mu }}_{\pm }|\mathrm{\Psi }).$$

Bayes' rule tells us the probability of a given observation Ψ given the probability of the parameters (marginal on the left or posterior) is proportional to the probability given the parameters (marginal on the right or prior) and the full distribution (likelihood),

$$\begin{array}{rl}& P(\mathit{\nu },\mathit{\lambda },{\mathit{\beta }}_{\mathrm{0}},\mathit{\delta }\mathit{\beta },{\mathit{\mu }}_{\pm }|\mathrm{\Psi })\propto P(\mathrm{\Psi }|\mathit{\nu },\mathit{\lambda },{\mathit{\beta }}_{\mathrm{0}},\mathit{\delta }\mathit{\beta },{\mathit{\mu }}_{\pm })\\ & \cdot P(\mathit{\nu },\mathit{\lambda },{\mathit{\beta }}_{\mathrm{0}},\mathit{\delta }\mathit{\beta },{\mathit{\mu }}_{\pm }).\end{array}$$

Sampling different values from the parameters' prior distributions will give corresponding values for the posterior distributions. A Bayesian sampler chooses successive values that tend towards greater likelihood of the model, given the observed trajectory, and will converge towards an optimal fit. Conceptually, this is what an MCMC (Markov chain–Monte-Carlo) optimiser does (Bolstad, 2010). A widely used sampling algorithm is the Metropolis algorithm (Hastings, 1970; Bernardo and Smith, 2009), which we also use here. This algorithm has been implemented in many software packages (e.g. Salvatier et al. (2016); Carpenter et al. (2017)).

The sampling process is time consuming because the evaluation of the potential (to calculate $P\left(\mathrm{\Psi }\right|\mathit{\nu },\mathit{\lambda },{\mathit{\beta }}_{\mathrm{0}},\mathit{\delta }\mathit{\beta },{\mathit{\mu }}_{\pm })$) requires numerical integration (using a quadrature method), which is costly to evaluate (the exponential family of distributions cannot, in general, be evaluated analytically).

3.1.1 Prior distributions

The prior distribution of a parameter represents all the information known about that parameter before confrontation with the observed values (Bolstad, 2010). With ν and λ transform the AMOC state variable (Ψ) with a shift (λ) and a scaling (ν). The shift λ cannot exceed the normalisation ν, giving an upper bound on λ. In addition, we note the lower limit of the lower branch, meaning λ must be larger than this minimum value. Similarly, the scaling ν cannot be larger than the maximum value of the AMOC on the upper branch. We expect the linear parameterisation of α and β introduced in the previous section to be 𝒪(1).

We are nonetheless still faced with infinite support on the coefficients of the expansion of the parameters (β₀, δβ). We therefore transform β₀ and δβ, with support $(-\mathrm{\infty },\mathrm{\infty })$, using the arctan function to map to $(-\mathit{\pi }/\mathrm{2},\mathit{\pi }/\mathrm{2})$. After such a transformation, we can sample from the flat prior distribution on that interval with most of the probability mass on “reasonable” values (i.e. 𝒪(1)). The following prior distributions are used:

$$\begin{array}{ll}{\displaystyle }\mathit{\nu }& {\displaystyle }=U(min(\text{AMOC}),max(\text{AMOC}\left)\right),\\ {\displaystyle }\mathit{\lambda }& {\displaystyle }=U(min(\text{AMOC}),\mathit{\nu }),\\ {\displaystyle }{\mathit{\mu }}_{+}& {\displaystyle }=U({\mathit{\mu }}_{\mathsf{\text{s+}}},{\mathit{\mu }}_{\mathsf{\text{up}}}),\\ {\displaystyle }{\mathit{\mu }}_{-}& {\displaystyle }=U({\mathit{\mu }}_{\mathsf{\text{dn}}},{\mathit{\mu }}_{\mathsf{\text{s-}}}),\\ {\displaystyle }\tan \left({\mathit{\beta }}_{\mathrm{0}}\right)& {\displaystyle }=U(-\mathit{\pi }/\mathrm{2},\mathit{\pi }/\mathrm{2}),\\ {\displaystyle }\tan \left(\mathit{\delta }\mathit{\beta }\right)& {\displaystyle }=U(-\mathit{\pi }/\mathrm{2},\mathit{\pi }/\mathrm{2}),\end{array}$$

where min(AMOC) and max(AMOC) are the minimum and maximum values in an observed collapse trajectory. U is the uniform distribution on indicated intervals. We stipulate the interval values of the collapse points μ_± as being bounded by where the trajectories merge (μ_𝙪𝙥 and μ_𝙙𝙣) and the inner values (μ_𝙨- and μ_𝙨+) observed in the trajectories (within which bimodality is demanded; see Fig. 6). ¹

3.2 Fitting EMIC collapse trajectories

An AMOC collapse was induced in models of intermediate complexity in Rahmstorf et al. (2005) by applying a freshwater forcing to the North Atlantic subtropical gyre region that reduced the salinity in the subpolar gyre to its north. Six of these models have a 3D ocean components; in Fig. 6 the trajectories of those collapses are reproduced (right column, the freshwater flux is labelled μ) together with their numerical derivatives (left columns in the panels). In Table 1 the models are listed. The forcing values of μ are known and are the same for each climate model. Each model was run to equilibrium for each forcing value; there is therefore no explicit time dependence in the hysteresis loops shown. Both the AMOC strength and the forcing value are given in units of Sv (=10⁶ m s⁻¹). Note that the bifurcation points (μ_±) must lie within the range where the trajectories appear bimodal.

Prange et al. (2003)Goosse et al. (2001)Edwards and Marsh (2005)Rahmstorf and Willebrand (1995)Weaver et al. (2001)

Table 1 Overview of models used. Each data point is independent of the others because each is the result of a quasi-steady-state run. The number of data points for each model was regridded onto a uniform freshwater forcing range consisting of 300 points. The summary of the type of model component and references are taken from Rahmstorf et al. (2005).

Table 2 Overview of models, the estimated standard deviation with the upper branch fitted to a linear function (note that the original trajectories had already been smoothed), the ranges of μ_±, the location of the present day in the models, and whether the present day value is in the unimodal regime (+) or not (−). All values are given in units of Sv.

The trajectories are from the numerical Earth System Models (EMICs) Rahmstorf et al. (2005, Fig. 2, bottom panel). The numerical derivatives show where the AMOC changes quickest as a response to the change in freshwater forcing. Each model has two peaks where the changes are largest, one for each change between stable branches. These peaks are located at the repeller in between the two attractors (the stable branches). At the repeller only unstable solutions exist and the AMOC is driven to a stable solution away from these states.

Figure 6 Absolute values of numerical derivatives (left) from the trajectories of AMOC strength as function of freshwater forcing to the right (taken from Rahmstorf et al., 2005, Fig. 2, bottom panel, reproduced with permission from the publisher: American Geophysical Union). In red the upper branch is shown, and in blue the lower branch is shown. The left column shows data from Bremen, ECBilt-CLIO, and C-GOLDSTEIN. The right column shows data from MOM hor, MOM iso, and UVic. Vertical solid lines mark μ=0 (blue) and μ=0.2 (red); vertical dashed lines mark the chosen boundary values for μ_±. All values are given in units of Sv.

If no other mechanisms apart from the salt advection are important, we expect the bifurcation points to lie beyond the observed transition points because a noise-induced transition pushes the AMOC into the off state sooner. Note that although the collapse points are expected to lie before these peaks, low levels of noise will obscure this effect. The dashed lines indicate the regions where we will search for the optimum values of μ_±. These differ from the fixed 0 and 0.2 values chosen by (Rahmstorf et al., 2005), who also shifted the trajectories to align on these values.

Before fitting, the upper and lower branches were extended to the left and right to fill the space of $-\mathrm{0.2}<\mathit{\mu }<\mathrm{0.4}$. A linear fit was used to produce additional values of the corresponding branches (at the same density of those points already present). All models then occupy the same freshwater forcing space. This is desirable because not all models have a lower branch that is fully sampled (specifically, UVic). The lower branch was extended with a negative rate of increase if the lower branch was moving upwards with increasing μ (MOM hor and MOM iso).

Our main goal is to model the transition from the on branch to the off branch, i.e. the upper right half of the hysteresis curve, and not so much the dynamics that govern the lower branch. In addition, because of this assumption, other dynamics govern the lower branch, and our simple model has to be extended to account for those dynamics. We ignore the data on the lower branch before the collapse point so that the fits would not be influenced by these points. We expect the remaining points of the trajectory to be dominated by the salt-advection mechanism.

We start by identifying some characteristic points in the trajectories in Table 2. The σ (variance of the process) of the models is not given in Rahmstorf et al. (2005) or elsewhere in the literature but was estimated as the deviation with a fitted function to the left most the top branch. Note that smoothing was already applied in Rahmstorf et al. (2005), lowering the variance of the trajectories. As we want to fit the collapse trajectory as given, we use the variance as evident from the data. In principle, σ could also be estimated as a parameter in the Bayesian optimisation, but that would unnecessarily enlarge the search space. Note that the “off-state” of the AMOC in these models is not 0, but ∼2 Sv of AMOC strength. If the salt-advection mechanism were the only operative effect, we expect this value to be ≤0. If a reverse advection cell emerges as the lower hysteresis branch, this value is negative.

In Fig. 7 fitted distributions are shown (also tabulated in Table 3). As best-fit parameters, we choose the mean values of the marginal posterior distributions. The dashed grey line marks the position of the unstable solution (repeller) in between the two attractor branches that separate the two basins of attraction.

Figure 7 Estimated distributions under changing μ: (a, b, c) Bremen, ECBilt-CLIO, and C-GOLDSTEIN; (d, e, f) MOM hor, MOM iso, and UVic. Vertical dashed lines mark the chosen boundary values for μ_±, with solid lines showing the fit values. The dashed grey line indicates the local minimum in the distribution (trench). The top ± symbols indicate the sign of the discriminant D for the fitted distribution (+ for unimodal, − for bimodal). Distribution spreads have been inflated with a factor ν∕2 to make them visible. All values have are given in units of Sv.

Table 3 Mean values and standard deviations of parameters corresponding to the fitted functions in Fig. 7. The root-mean-square deviation (a goodness of fit measure) has been determined on the upper branch up to the fitted collapse point.

The fits with a linear series through the (α,β) parameter space result in a mismatch between the behaviour seen on lower branches and that on the upper branches. This is less obvious for UVic and ECBilt-CLIO but is especially apparent for the two MOM models.

4 Discussion and conclusion

We derived a simple model of AMOC collapse based on Langevin dynamics (Eq. 1) with a changing freshwater forcing (μ) and applied this to EMIC-simulated collapse trajectories taken from Rahmstorf et al. (2005). The collapse occurs at a bifurcation point μ₊ that appears smaller than given in (Rahmstorf et al., 2005). A corresponding bifurcation point μ₋ relates an abrupt transition back to the on state. The AMOC also requires an offset and scaling parameter to be fitted (λ and ν). These six parameters are sufficient to describe the abrupt collapse of the AMOC as part of a hysteresis loop under varying freshwater forcing.

Any process that allows two stable states with rapid transitions between them and an asymmetric response to the forcing could in principle be described by our method. Other such geophysical processes can include ice sheet mass loss (e.g. Robinson et al., 2012), forest dieback (e.g. Staal et al., 2016), or lake turbidity (Scheffer and van Nes, 2007).

The resurgences of the AMOC seen in the hysteresis diagrams behave differently from the collapses. The Langevin model is too simple to capture both processes. It is, however, possible to fit the change in the upper branch of the AMOC – the on state – as it moves towards a critical point and the dominant salt-advection feedback mechanism breaks down.

We note that Rahmstorf et al. (2005) determine the AMOC strength as the maximum of the meridional volume transport in the North Atlantic, and this might explain the asymmetry between the two branches. If for a reverse overturning cell the wrong metric has been used, then the lower branch location is not correct. It is conceivable that the Langevin model results in better fits if Rahmstorf et al. (2005) had sampled $max\left(\right|\mathrm{\Psi }\left|\right)$ instead of max(Ψ), which would have resulted in a better metric of the lower branch. With the metric used it is not apparent whether a reversed overturning cell was present or not because it was not sampled if the AMOC had taken on a negative value. It is unclear to what extent the models discussed here develop a reversed overturning circulation, which can arise in 3D models (Weijer and Dijkstra, 2001; Yin and Stouffer, 2007) but can also be suppressed by atmospheric feedbacks (Yin and Stouffer, 2007; however, see also Mecking et al., 2016) and strongly affected by gyre dynamics (Prange et al., 2003). These effects are not captured by the simple Langevin model proposed here, but at present it is still unclear to what extent these effects are essential in capturing the first-order stability properties of the AMOC. In each case, there is no obvious way to model the asymmetry between the two branches and obtain a full description. The two branches could be separated by associating each with a different overturning cell. The upper branch is identified with the NADW-driven cell, while a reverse cell is responsible for the lower branch. If a reverse overturning cell (as described in e.g. Yin and Stouffer, 2007) indeed dominates the lower AMOC branch, two separate overturning cells are responsible for the observed trajectories, and the two branches then cannot be expected to fit with the same parameter set.

However, another possible explanation is that (two) separate mechanisms are responsible for the upper and lower branch dependency on μ. Possible mechanisms include possible mechanisms include the influence of wind stress, North Atlantic subpolar gyre convective instability (Hofmann and Rahmstorf, 2009), or other pathways of deep water formation (Heuzé, 2017). In addition, changes in the ITCZ (inter-tropical convergence zone) due to ocean-atmosphere feedbacks are possible (Green et al., 2019); these can, in turn, can affect the salinity of the North Atlantic subtropical gyre region. However, Mecking et al. (2017) showed that for a high-resolution model the salt-advection feedback was nevertheless stronger than the ITCZ effects. Other wind coupling can occur further south through a coupling with the ACC (Antarctic Circumpolar Current) which is based on the thermal wind relation (Marshall and Johnson, 2017).

A third explanation is that deep water formation is a local process, and as a result an asymmetry is to be expected between the two branches. Local convection can, however, be subject to global controls and be associated with a sinking branch which occurs in conjunction with deep convection, but is not directly driven by it, see Spall and Pickart (2001) for a detailed discussion. The AMOC could develop a reverse cell where the overturning is driven by Antarctic Intermediate Water (AAIW), which is not part of the conceptual picture presented here (Yin and Stouffer, 2007; Jackson et al., 2017). The reverse cell introduces an asymmetry in the collapse trajectories because the driver of deep water formation is not in the North Atlantic, and might break our assumption that both the on and off branches are controlled by the same process. It is therefore difficult to estimate the return path of the AMOC if the lower branch has additional drivers from the dominant salt-advection mechanism of the upper branch. Forcing values appropriate for the lower branch might be different than those found for the upper branch.

Furthermore, the methodology used in this paper comes with difficulties in the numerical implementation. The fit procedure requires the normalisation of each distribution in the μ time series. Because no analytic solution exist a numerical approach is needed. The numerical integration adds to the computational costs of the fits. The Markov chain method is also prone to find local optima. In addition, the cost of numerical integration necessitates stopping the fits at shorter chains than (perhaps) are needed, an analytic formulation of the integrand would alleviate this but none exists to our knowledge. Modern sampling algorithms allow for gradient information to be used, which is effective when sampling a higher dimensional parameter space (the Metropolis algorithm used in this paper has greater difficulty as the dimensionality of the parameter space increases). Tighter constraints on the prior distributions could be beneficial here.

As stated in Rahmstorf et al. (2005), the EMIC trajectories had already been smoothed, resulting in a smaller variance; a smaller variance leads to distributions that are more sharply peaked. This increases the computational cost of integrating the distributions numerically. Smoothing can also add to the inertia seen in the collapses, but might be due to other reasons such as stopping the EMIC simulations before equilibration of the AMOC collapse, leaving the AMOC in a winding-down state. In addition, the models in Rahmstorf et al. (2005) were integrated for 1000 model years per freshwater forcing value (which was changed in 0.05 Sv increments). If the integrations were done for an insufficient amount of time, the AMOC collapse is incomplete, leaving the measured value out of equilibrium. The intermediate points in the collapse trajectories beyond the bifurcation points indicate that either the sample points are inaccurate or other processes are involved in the AMOC.

Finally, the fitted collapse trajectories were done on an ensemble of EMICs, which arguably are not sufficiently representative of the real climate. As noted by Gent (2018), the hysteresis behaviour has not been investigated fully in models of greater complexity than EMICs; the computational cost being prohibitive for models with high resolution (and short time steps). The hysteresis behaviour in glacial state changes has, however, been investigated in greater detail using models with simplified dynamics (e.g. Schiller et al., 1997; Zhang et al., 2017). The question arises to what extent the procedure outlined in this paper can be applied to more complicated models such as those in the CMIP archives (Taylor et al., 2012). These models do not show a full collapse trajectory like those in Rahmstorf et al. (2005), which means no sample points of the lower branch are available. Also, CMIP provides times series of forced runs. To validate our method, a transient run requires known equilibrium bifurcation points under a slowly changing μ and include an AMOC collapse. Using a simple box model, transition probabilities for an AMOC collapse have been determined by Castellana et al. (2019). From the CMIP ensemble a similar estimate might be obtained, or at least the collapse characteristics of various models can be compared. Provided the CMIP models accurately capture the behaviour of the real AMOC and the freshwater forcing counterpart (our μ) can be identified, an estimate can be made of the distance of the current climate state to the collapse point. Freshwater quantities such as M_ov have been posited (e.g. Drijfhout et al., 2011) as being suitable indicators of AMOC stability. It is possible that M_ov relates to μ and can be used to extend our method to transient runs, but at present it is unknown whether this can be done. The inclusion of ice sheets can make a substantial difference in AMOC recovery (Ackermann et al., 2020). In addition, the atmospheric freshwater transport might have a stabilising effect on the AMOC that is greater than the freshwater transports by the ocean (Lohmann, 2003). There is, however, also evidence that coupled climate models suffer from a salinity bias that favours an AMOC that is too stable (Drijfhout et al., 2011; Liu et al., 2017). These matters are outside the conceptual picture of M_ov as a stability indicator. It is therefore still an open question how probable an AMOC collapse is in more realistic models and in reality. However, using the method outlined in this paper, a first step could be made towards answering this question.