Collapse of the Atlantic Meridional Overturning described by Langevin dynamics

Using a machine learning technique, collapse trajectories of the Atlantic Meridional Overturning Circulation 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 under a freshwater forcing. Reversing the freshwater forcing results in asymmetric behaviour that is less well captured and would require a more complicated model.

minimum required for having three distinct solutions (double wells). This potential function has two parameters which arepresumed 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.
Though the Langevin equation has played a role in the conceptual picture of bistability and tipping points in the climate, but it has not been used to actually fit the parameters to a (simulated) AMOC collapse. Here, we attempt to construct a simple 30 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.
Section 2 sketches the theoretical background of the Langevin equation and of the salt-advection mechanism. In Section 3 we fit the proposed Langevin model to the AMOC collapse trajectories seen in a set of climate models of intermediate complexity 35 (EMICs) taken from Rahmstorf et al. (2005). We end with a discussion and conclusions in Section 4.

The Langevin model
An increase in surface air temperatures, or an increase in P-E (precipitation -evaporation) of freshwater surface flux, 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 reduced. In principle, this mechanism can reduce the AMOC 40 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, the deep southward flow couples to the surface return flow via upwelling or other mechanisms in the Southern Ocean, affecting the salinity in the North Atlantic subpolar gyre. A reduction in salinity decreases buoyancy, and this positive feedback accelerates the process of AMOC weakening and a collapse results on relatively short timescales. 45 Fig. 1 shows a conceptual picture of the two stable AMOC (index) states. The AMOC is a zero-dimensional variable arrived at 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 50 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 ceases to exist and a qualitative change takes place in the potential function (the number of solutions for a given value of the freshwater forcing µ goes from 2 to 3).
Below we will derive a model based on the Langevin equation that captures the essential dynamics of a bimodal AMOC under a freshwater forcing µ.  3 https://doi.org/10.5194/esd-2020-43 Preprint. Discussion started: 30 June 2020 c Author(s) 2020. CC BY 4.0 License.

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), 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 (deterministic) bistability seen in Fig. 1 has been studied mainly in a qualitative way (within catastrophe theory). In Poston and Stewart (1978) an extensive treatment is given why only two parameters are sufficient to describe the bistability.

65
These two parameters describe the critical behaviour, not just locally near the critical points, but the entire trajectory under a suitable transformation. (The behaviour at small scale is fundamentally tied to the global behaviour.) 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 form (Gardiner, 2004;Ditlevsen and Johnsen, 2010) The two parameter 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 ν. The value of ν is roughly the distance in Ψ between the bifurcation point on the top branch to the bifurcation point on the lower branch. The value of λ is roughly 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 80 the demand that the freshwater forcing is the only variable that determines the dynamical behaviour.

Potential description
In Fig. 2  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 landscape is shown where the areas indicated are those with qualitatively different behaviour seen in Fig. 2. See also Poston and Stewart (1978) for similar diagrammes. The cusp point P is the singular point where no proper 90 solution can exist because only the trivial solution is allowed here (both collapse 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 (µ).
Our aim is to arrive at a description that matches a track of µ values across the stability landscape.The two parameters α, β are independent but can be parametrised by other variables that map them to observations. In the literature α is referred to as 95 the normal factor, and β the splitting factor (Poston and Stewart, 1978). If parametrised by a single variable the track across the stability surface is one-dimensional, 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.

Constraints
With a varying α there exist two critical points (α ± ) in between which the distribution is bimodal and unimodal outside that 100 interval. Because the AMOC trajectory is 1-dimensional, there must be a 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 (real solutions) needs to be solved for α to obtain the two critical solutions 105 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 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 α ± and β ± the values corresponding to µ ± . The track across the stability landscape is 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 α ± there is a corresponding β ± at the respective critical points) and a varying β 120 corresponds to a slanted track across the stability landscape in Fig. 3. We assume linear functions for α and β, reducing the dependency to these four parameters. Linear functions are the most simple non-trivial dependencies, while adding non-linear parameters introduces further unknowns, making this the most parsimonious parametrisation that captures the first 125 order behaviour. Also, intuitively we can understand the pair (α, β) as the angle under which the system moves to the bifurcation 7 https://doi.org/10.5194/esd-2020-43 Preprint. Discussion started: 30 June 2020 c Author(s) 2020. CC BY 4.0 License. point (B 1,2 ) in Fig. 3), which locally only requires the values of α and β up to first order. Poston and Stewart (1978, p. 59) also remark that the system's local behaviour is essentially the same between critical points, which means a linear expansion should suffice for fitting the upper branch. 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 parametrised by β 0 and δβ. Note that only solutions with β ± > 0 are valid. Also, values for β 0 and δβ 0 that result in a track that crosses B 2 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.)

140
With the deterministic framework in place, the stochastic nature can be reintroduced.We obtain a distribution needed to fit the parameters of the potential function. As shown by Cobb (1978), this distribution belongs to the exponential family.
For a polynomial function as the potential, the distribution obtained is from the exponential family, which is a generalisation of the exponential distribution (Balakrishnan and Nevzorov, 2004) where any function can determine the exponent value. The potential, we had already supposed to be described by a fourth-order polynomial in the previous section, gives the probability Note that C = C(α, β) and does not have a (known) analytical expression for the general case. The value of σ is a measure of intrinsic variation in the AMOC. See Gardiner (2004) for a derivation of this distribution using the Fokker-Planck equation, from which also the Langevin equation can be derived. Note that σ→σ/ν because to the scaling with ν we introduced in the AMOC rapidly moves from one attractor to the other. For this reason, the sampled bimodality region might be larger than is apparent from a particular sample AMOC trajectory.
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. These changes correspond directly to the potential functions 160 in Fig. 2. In principle, a change in on only one of the two structural parameters (α and β) can move the distribution between unimodal and bimodal forms.
The normalisation of the family of distributions depends on the values of the parameters. Therefore, we are required to calculate the normalisations factors for each parameter set. This cannot be done analytically, but can be done accurately with an adaptive quadrature method (Piessens et al., 2012), though it suffers from numerical limitations. in the middle and inversion from weak to dominant takes place. The relative strength between dominant and weak depends on σ. Right: 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.
We are now in a position to apply the above to collapse trajectories from climate models.

Simulated AMOC collapse parameter estimation
We describe how to find an optimal solution under the framework arrived at 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 estimate the values for µ ± , resulting in a six parameter list that describes (the upper branch) of an AMOC collapse.

170
The parameters β 0 and δβ are independent to each other, but need to cross the curves B 1,2 in Fig. 3  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 https://doi.org/10.5194/esd-2020-43 Preprint. Discussion started: 30 June 2020 c Author(s) 2020. CC BY 4.0 License.
trajectory parameters β 0 and δβ requires a more flexible method. Knowing which distribution to use, we can fit its parameters under some measure of goodness-of-fit by machine learning. Specifically, we can use Bayes' rule to maximise the likelihood 180 of a parametrised Langevin model L σ (ν, λ, β 0 , δβ, µ ± ) given a trajectory Ψ(µ) (Bolstad, 2010), to arrive at the (linearised) posterior distributions of ν, λ, β 0 , δβ, µ ± under the observed Ψ(µ). 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), does (Bolstad, 2010). A widely used sampling algorithm is the Metropolis algorithm (Hastings, 1970), which we also use here. 1 The models can be fit with uninformative priors, but the sampling process is time consuming because the evaluation of the potential (to calculate P (Ψ | ν, λ, β 0 , δβ, µ ± )) requires numerical integration (using a quadrature method), which is costly to 195 evaluate (the exponential family of distributions cannot, in general, be evaluated analytically).

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 λ introduced earlier, the state variable x undergoes an affine transformation and normalises the polynomial. These transform the AMOC state variable (Ψ) with a shift (λ) and a scaling (ν). The shift λ cannot 200 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 AMOC on the upper branch. We expect the linear parametrisation of α and β introduced in the previous section to be O(1).

Fitting EMIC collapse trajectories
An AMOC collapse was induced in six models of intermediate complexity in Rahmstorf et al. (2005) by applying a freshwater 220 forcing to the North-Atlantic subtropical gyre region that reduced the salinity in the subpolar gyre to its north. In Fig. 6 Table 1. 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 present day in the models, and whether the present day value is in the unimodal regime (+) or not (-). All values in units of Sv.
The trajectories are taken from the numerical Earth System Models (EMICs) Rahmstorf et al. (2005, Fig. 2, bottom panel) 225 by extracting the data points directly from the graphic in the electronic publication. 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 repellor in between the two attractors (the stable branches). At the repellor only unstable solutions exist and the AMOC is driven to a more stable solution, away from these states. 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. The dragged-out descent to the lower branch (e.g. the model 235 Bremen) indicates that the salt-advection mechanism does not necessarily result in an abrupt collapse in the trajectory.
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 240 moving upwards with increasing µ (MOM hor and MOM iso).
Because we want to model the collapse of the AMOC, only the behaviour of the upper branch and the lower branch beyond µ + is of interest. We expect these points of the trajectory do be dominated by the salt-advection mechanism. The hysteresis loop is exemplified in the trajectory of C-GOLDSTEIN. This model shows a smooth path with a relatively rapid transition region-but not a collapse as such by appearance. 245 We start by identifying some characteristic points in the trajectories in Tab. 1. 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 (part of) the top branch (note that smoothing was already applied in Rahmstorf et al. (2005), lowering the variance of the trajectories). Note that the 'off-state' of the AMOC in these models is not 0, but ∼ 2Sv 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 250 hysteresis branch, this value is negative.
The model trajectories are apparently driven by the forcing value, which means we should be able to explain the behaviour of the hysteresis using this variable as the only driver applied to the Langevin model.
In Fig. 7 fitted distributions are shown for linearly parametrised sample paths through the stability space (also tabulated in Tab. 7). The β parameter changes linearly and α follows from the constraints in Eqs. 4 and 5. Blue and red lines indicate 255 the prior bounds for µ − and µ + , respectively. The dashed grey line marks the positions of the unstable solution (repellor) in between the two attractor branches which separates the two basins of attraction.
The fits with a linear track through the (α, β) parameter space result in a mismatch between the behaviour seen on lower branches and that on the upper branches. In particular, the upper branch shows a non-linear degradation, where the fitted distributions do not. This is less obvious for UVic and ECBilt-CLIO, but especially apparent for the two MOM models.

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 µ + which appears smaller than given in (Rahmstorf et al., 2005). A corresponding bifurcation point µ − relates an abrupt transition back to the on-state. Additionally, a linear parameterisation through state space couples to the freshwater applied to 265 the North Atlantic subtropical gyre region. The AMOC also requires an offset and scaling parameter to be fitted (λ and ν).
These six parameters are sufficient to describe the abrupt collapse/re-invigouration of the AMOC that leads to a hysteresis loop under varying freshwater forcing.
The AMOC collapse and re-invigouration seen in these models cannot be completely fitted with Langevin dynamics due to the asymmetry in the lower vs the upper branch. It is, however, possible to fit the change in the upper branch of the AMOC-the We note that Rahmstorf et al. (2005) determine the AMOC strength as the maximum of the meridional volume transport in the North Atlantic and 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(Ψ)| instead of max(Ψ), which would have resulted in a better metric of the 275 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. Therefore, 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 indeed a reverse overturning cell (as described in e.g. Yin and Stouffer (2007)) dominates the lower AMOC branch, 280 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, SPG convective instability (Hofmann and Rahmstorf, 2009), or other pathways of deep water formation (Heuzé, 2017). Also, changes in the 285 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. Also, Gent (2018) notes that the EMICs in Rahmstorf et al. (2005) have reduced air-sea interaction feedbacks compared to more modern and more complicated fully coupled ocean-atmosphere models; stronger feedbacks lead to greater AMOC stability. Other wind coupling 290 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 detailed discussion.

295
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 300 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 needed. The numerical integration adds to the computational costs of the fits. The Markov chain method is prone to find local optima. Also, the cost of numerical integration necessitates stopping the fits at shorter chains than (perhaps) are needed, an 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). It is possible for non-admissible solutions to be generated; in which case, the sampler is effectively reduced to a random searcher, till its finds a solution subspace to optimise under. (Tighter constraints on the prior distributions could be beneficial here.)

310
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. Also, it is not clear how long the models in Rahmstorf et al. (2005) were integrated per freshwater forcing value. If the 315 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. The question arises to what extent the procedure outlined in this paper can be applied to more complicated 320 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. If it is indeed possible to use direct numerical stochastic integration of the Langevin equation, no lower branch needs to be sampled. Sampling from the prior distributions and optimising under the observed upper branch should also lead to robust estimates. Another advantage is that the probability distribution can be assumed Gaussian when on the upper branch. Using a simple box model, transition 325 probabilities for an AMOC collapse have been determined by Castellana et al. (2019). From the CMIP ensemble a similar estimate might be obtained using a direct numerical stochastic integration approach. 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. It is therefore still an open question how probable an AMOC collapse is in more realistic models, and reality, but with the method outlined in this paper a first step could be made Competing interests. The authors declare that they have no conflict of interest.