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.

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 (repellor) 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 60 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 µ.

Multiple stable AMOC states
The conceptual picture of the AMOC being a zero-dimensional variable that is driven by stochastic forces trapped in a po-65 tential 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), 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 70 and focus on the deterministic behaviour.
The double well potential seen in Fig. 1 has been extensively studied and applied, also in a quantitative way. But to our knowledge it has not been quantitatively applied to AMOC hysteresis using the Langevin equation in complex numerical climate models before.
studied mainly qualitatively in connection with the Langevin equation. AMOC bistability has, however, been studied quan-75 titatively 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. From a fourth order polynomial for U ::::: More :::::::: precisely, the third and fourth order coefficients :::: 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 the entire trajectory under a suitable transformation. A direct consequence is that only partial information, in the form of a piece 80 of the trajectory, should suffice to describe the entire trajectory (the full hysteresis loop).
The two parameters α, β are functions of the freshwater forcing µ. The AMOC state variable Ψ requires an affine transformation (Cobb, 1980), 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 β the splitting factor (Poston and Stewart, 1978). In the 90 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. Below these are the three possible unimodal states (E). These occur for forcing values to the left of µ − and to the right of µ + .

100
The panels B 1 and B 2 are the submanifolds 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.
In Fig. 3 the stability diagram is shown where the areas indicated are those with qualitatively different behaviour seen in Fig. 2. See also 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 105 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 (µ).
Our aim is to arrive at a description that matches a series of µ values across the stability diagram. The two parameters α, β 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 110 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.

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 1-dimensional and µ is also 1-dimensional, there must be a 115 relation between α and β that reduces dimensionality from two to one dimensions. 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 ∂ x U , µ ± , can be found analytically for µ ± real and being degenerate solutions.
It can be shown (Birkhoff and Mac Lane, 1970, p. 106) that the discriminant D = 27α 2 − 4β 3 = 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 120 changes at forcing values µ ± . (When D < 0 there are three distinct real solutions which corresponds 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 3 to 2 (or 1 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 β, 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 α + < 0 −U (1) < 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 with α ± and β ± the values corresponding to µ ± . Changing µ in the bifurcation diagram corresponds to moving from curve B 2 to curve B 1 and Eq. 3 relates the two stability parameters α and β at the two critical forcing values µ ± .

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 135 assume linear functions for α and β, 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 parametrisation that captures the first 140 order behaviour. Also, 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 parametrisation we can determine the offset α 0 and rate δα in terms of β 0 and δβ, This constrains the values of α, leaving only β as a free variable, which is then parameterised by β 0 and δβ. Note that only solutions with β ± > 0 are valid. Also, values for β 0 and δβ 0 that result in crossing B 2 in another point besides β − are 150 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.)

Stochastic interpretation
With the deterministic framework in place, the stochastic nature can be reintroduced.The potential function can be replaced by a distribution which is the stationary distribution in the asymptotic limit (i.e. the long term behaviour of repeated sampling of 155 the hysteresis loop). The potential (a fourth-order polynomial) gives the probability distribution (Cobb, 1978) P (x, α, β) = Ce −2/σ 2 U (x) = Ce 2/σ 2 (−1/4x 4 +β/2x 2 +αx) .
The factor C = C(α, β) 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 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 also the Langevin equation can be derived. Also, note that σ→σ/ν because of the scaling with ν we introduced in Section 2.1. The distributions in Figs 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 175 the distribution between unimodal and bimodal forms.
We are now in a position to apply the above to collapse trajectories from climate models.

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 β 0 and δβ can be found, together with the scaling parameters ν and λ. We will also 180 estimate the values for µ ± , resulting in a six parameter list that describes (the upper branch) of an AMOC collapse.
The parameters β 0 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 α 0 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 185 are located already gives a rough estimate. 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 β 0 and δβ requires a more flexible method. Knowing which distribution to use, we can estimate the 190 posterior probability distribution of the parameters given the data Ψ(µ),

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 λ. Also, 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 210 AMOC on the upper branch. We expect the linear parametrisation of α and β introduced in the previous section to be O(1).
We are nonetheless still faced with infinite support on the coefficients of the expansion of the parameters (β 0 , δβ). We therefore transform β 0 and δβ, with support (−∞, ∞), using the arctan function to map to (−π/2, π/2). After such a transfor-mation, we can sample from the flat prior distribution on that interval with most of the probability mass on 'reasonable' values (i.e. O(1)). The following prior distributions are used: ECBilt-CLIO 243 3D primitive equations quasi-geostrophic Goosse et al. (2001) C-GOLDSTEIN 849 3D simplified energy-moisture balance Edwards and Marsh (2005) MOM hor 1233 3D primitive equations (MOM) simple energy balance Rahmstorf and Willebrand (1995) MOM iso 1442 as above, with isopycnal mixing simple energy balance UVic 464 3D primitive equations (MOM) energy-moisture balance Weaver et al. (2001) Table 1. Overview of models used. Each data point is independent from the others because each is the result of a quasi steady state run.
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.

245
Before fitting, the upper and lower branches were extended to the left and right to fill the space of −0.2 < µ < 0.4. A linear fit was use 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).

250
Our main goal is to model the transition from on-branch to :: the : off-branch, that is, the upper right half of the hysteresis curve, and not so much the dynamics that govern the lower branch. Also, because we assume that 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 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. 255 We start by identifying some characteristic points in the trajectories in Tab. 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. Because 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.
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.

Discussion and conclusion 270
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 µ + which 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 that leads to :: as :::: part :: of : a hysteresis loop under varying 275 freshwater forcing. The resurgence of the AMOC is not the same as the collapse process and we did not attempt to obtain an accurate fit of that part of the hysteresis loop.
Any process which 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 might be ice sheet mass loss (e.g. Robinson et al. (2012)), forest dieback (e.g. Staal et al. (2016)), and lake turbidity (Scheffer and van Nes, 2007).
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). Also, changes in the ITCZ (inter-tropical convergence zone) due to ocean-atmosphere feedbacks are possible (Green 305 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 310 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 315 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 µ timeseries. Because no analytic solution exist a numerical approach is 320 needed. The numerical integration adds to the computational costs of the fits. The Markov chain method is also prone to find local optima. Also, 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 325 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.
Comments to the Author: Dear authors, the manuscript has improved and requires only minor revisions.
1) please follow the suggestions of the referee, and provide a point-by-point answer to the comments.
We thank the editor for the additional comments and suggestions for improvement. Below our responses.
2) My own evaluation and following comments can improve the quality of the paper. a) page 2, lines 25-28: Referencing to double well potential and dynamics of abrupt climate shifts. I think the method and application of Livina and others were first described here: G. Lohmann, 2009: Phys. Rev. E, doi: 10.1103/PhysRevE.80.066104 This paper predates the reference to Livina et al (2010) and relates to ice core data. A stochastically driven motion in a 4 th order polynomial potential is used to model glacial state changes on the millennial scale. We have included it as an additional reference.
b) I think that the statement "1 This algorithm has been implemented in many software packages." could be avoided as footnote. You could include this statement in the full text and provide references.
We have included references to two implementations in the main body of the text. c) section4, lines 270ff, 282 ff. The Langevin model is a simplified view and many features like the gyres are missing. In Prange et al. (2003) it has been shown that even 2D and 3D models can have a substantial different dynamics because of the missing dimension. I guess that this is even more true for the 0-d model. In 3 D models, a reversed flow as off mode does not exist (see e.g. on Fig. 9 of their paper) due to the gyres (section 4b). THis might be also related to your statements lines 250 ff on page 13, and lines 260 ff on page 14, lines 274 and 283 on page 275. The salinity advection mechanism is indeed not the mechanism for the off state. If feel that the speculation with the metric (lines 284-286) could take this into account. The way it is written now is somehow misleading. Indeed, the statements on lines 250 and 260 relate to the differences in upper and lower branches. A lower dimensional model necessarily needs to ignore geometry that can be found in higher dimensional ones. We have added a remark to highlight this. "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 which 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 effect 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..." d) page 17, line 322: This statement is not correct: "it is not clear how long the models in Rahmstorf et al. (2005) were integrated per freshwater forcing value". This value is given in the paper. "The rate of change of the freshwater input was 0.05 Sv per 1,000 model years." This information could be easily included.
We added the integration times and forcing value increment value. e) page 18, lines 327 ff: The statement "As noted by Gent (2018), the hysteresis behaviour has not been investigated fully in models of greater complexity than EMICs" is probably not true. See some