The Effect of Indian Ocean Surface Freshwater Flux Biases On the Multi-Stable Regime of the AMOC

Original Research Papers

Authors

Henk A. Dijkstra
René M. van Westen

Abstract

It is known that global climate models (GCMs) have substantial biases in the surface freshwater flux. Using numerical bifurcation analyses on a global ocean model, we study here the effect of freshwater flux biases over the Indian Ocean on the multiple equilibrium regime of the Atlantic Meridional Overturning Circulation (AMOC). We find that a (positive) freshwater flux bias over the Indian Ocean shifts the multiple equilibrium regime to larger values of North Atlantic freshwater input but hardly affects the associated hysteresis width. The magnitude of this shift depends on the way the anomalous North Atlantic freshwater flux is compensated. The changes in bifurcation diagrams can be explained from the overall freshwater balance over the Atlantic basin.

Keywords:

Year: 2024

Volume: 76 Issue: 1

Page/Article: 90–100

DOI: 10.16993/tellusa.3246

Submitted on Jul 13, 2023

Accepted on Apr 17, 2024

Published on May 6, 2024

Peer Reviewed

CC Attribution 4.0

1 Introduction

The Atlantic Meridional Overturning Circulation (AMOC) has been proposed as one of the tipping elements in the climate system (Lenton et al., 2008; Armstrong McKay et al., 2022), indicating that it may undergo a relatively rapid change under a slowly developing forcing. The AMOC is thought to be particularly sensitive to the freshwater forcing, either through the surface freshwater flux, through input of freshwater due to ice melt (e.g. from the Greenland Ice Sheet) or through river runoff (Rahmstorf et al., 2005). These freshwater fluxes affect the meridional density differences in the ocean which control the AMOC strength and change the heat and salt transport carried by AMOC (Marotzke, 2000).

From conceptual models, the tipping behaviour is clearly related to the multi-stability properties of the AMOC. For example, in the Stommel two-box model (Stommel, 1961), there is an interval of the surface freshwater forcing where two stable steady AMOC states exist and tipping occurs due to transitions between these states. The central feedback responsible for the tipping behaviour is the salt-advection feedback, where a freshwater perturbation in the North Atlantic causes a weakening of the AMOC which leads to less northward salt transport and hence amplification of the perturbation (Marotzke, 2000).

Precise boundaries of the multiple-equilibrium regime of the AMOC have been obtained using conceptual models (Cessi, 1994; Cimatoribus et al., 2012) and fully-implicit ocean-climate models (De Niet et al., 2007; Toom et al., 2012; Mulder et al., 2021). One of the important results of these studies is that the existence of the multiple-equilibrium regime can be related (in these models) to an observable quantity (Rahmstorf, 1996), now often called the AMOC stability (or regime) indicator. This indicator has many different notations in the literature, e.g. M_ov (De Vries and Weber, 2005) or F_ov (Hawkins et al., 2011). Here, we will follow Weijer et al. (2019) and use F_ovS (F_ovN) as the freshwater transport carried by the AMOC over the southern (northern) boundary at 35°S (60°N) of the Atlantic basin (Dijkstra, 2007; Huisman et al., 2010; Liu et al., 2017). Available observations (Bryden et al., 2011) show that the present-day AMOC is exporting freshwater out of the Atlantic (F_ovS < 0). It is known that F_ovS ignores some relevant processes (Gent, 2018), but if one accepts that F_ovS is a proper indicator, the AMOC is in a multiple-equilibrium regime based on its observed values (Weijer et al., 2019).

Less precise estimates of the multiple-equilibrium regime boundaries can be obtained from global climate models. In so-called quasi-equilibrium experiments, the freshwater forcing is changed very slowly such that the model state stays close to the (slowly changing) equilibrium. When the freshwater forcing is varied in both directions and covers the multiple-equilibrium regime, regime boundaries can be inferred from the so-called hysteresis width, i.e. the freshwater forcing values where the AMOC collapses and recovers. The rate of forcing is important here and if this is much faster than the equilibration time scale of the steady state, the approximations of the regime boundaries become worse and also rate-induced tipping may occur (Lohmann et al., 2021). Such quasi-equilibrium experiments have been performed with many ocean-only models (Cini et al., 2024; Lohmann et al., 2024), Earth System Models of Intermediate Complexity (EMICs) (Rahmstorf et al., 2005) and the FAMOUS model, the latter being a Global Climate Model (GCM) with a relatively coarse horizontal ocean resolution of 2.5° × 3.75° (Hawkins et al., 2011). AMOC hysteresis behaviour has also been investigated in the Community Climate System Model (CCSM3) (Hu et al., 2012) and recently in the Community Earth System Model (van Westen and Dijkstra, 2023; van Westen et al., 2024).

Mostly due to computational constraints, the AMOC response to only particular freshwater forcing perturbations is often considered in state-of-the-art GCMs. In these so-called ‘hosing experiments’ (Stouffer et al., 2006), quite a diversity of model behaviour is found. It is not known whether a multiple-equilibrium AMOC regime exists in such models and only sporadic indications of such a regime have been found (Mecking et al., 2016; Jackson and Wood, 2018a, b). The problem is that it is difficult to assess whether the weak AMOC states resulting from the freshwater input are equilibrium solutions of the models. What these model studies certainly have shown is that an AMOC weakening would have severe impacts on the climate system, affecting sea level and regional temperatures in many areas around the world (Vellinga et al., 2002; Jackson et al., 2015; Liu et al., 2020; Orihuela-Pinto et al., 2022).

Of course, the real present-day AMOC may not have a multiple equilibrium regime and GCMs may model that correctly. On the other hand, the real AMOC may be in a multiple equilibrium regime, and the GCMs may not capture it. In that case, the GCMs misrepresent (or miss) crucial processes, such that they do not display tipping behaviour. A prominent example is the incorrect representation of ocean eddy transport processes in GCMs, which may prevent the existence of a multi-stable AMOC regime. Another possibility is that GCMs capture the relevant processes but the parameters in the models are not correct, such that a multi-stable regime does not occur due to model biases. An example of this is that many GCMs are considered to have a too stable AMOC due to biases in the freshwater transport in the Atlantic Ocean (Drijfhout et al., 2013; Mecking et al., 2016).

In van Westen and Dijkstra (2024), surface salinity biases were evaluated in the Coupled Model Intercomparison Projects (CMIP) phase 6 (CMIP6) models using historical simulations and comparing them to reanalysis data over the period 1994–2020. It was shown that several persistent biases in state-of-the-art climate models lead to an AMOC with an Atlantic freshwater transport that is in disagreement with observations (i.e., many models have F_ovS > 0). They also demonstrate that increasing the horizontal ocean and atmospheric resolution in the Community Earth System Model (CESM) does not reduce the various F_ovS biases. The largest F_ovS bias is caused by a too positive surface freshwater flux over the Indian Ocean, which is shown in Figure 1 for the high-resolution CESM (HR-CESM, 0.1° ocean and 0.25° atmosphere) and low-resolution CESM (LR-CESM, 1° ocean and 1° atmosphere) compared to reanalysis (i.e., ERA5). These positive freshwater fluxes cause a too fresh Indian Ocean and (via advection) Agulhas Leakage. This leads to an overestimation of the Atlantic freshwater import at 34°S and hence to larger values of F_ovS than in reanalysis.

Figure 1

The present-day (1994–2020) P-E bias (w.r.t. ERA5) for the (a): HR-CESM and (b): LR-CESM. The spatially-averaged P-E biases over the red-outlined region are 0.04 Sv and 0.05 Sv, respectively.

Motivated by these results, we extend the bifurcation analysis on the global ocean-climate model as in Dijkstra (2007) and determine the effects of Indian Ocean surface freshwater flux biases on the multi-stable regime of the AMOC. In section 2, the fully-implicit model used is shortly summarised and the continuation methodology to compute bifurcation diagrams is described. Then in section 3, we focus on the effect of freshwater biases on the bifurcation diagrams. Mechanisms of the shift in the multi-stable regimes are analysed using the freshwater balance over the Atlantic. A summary and discussion follows in section 4.

2 Formulation

2.1 Model

The fully-implicit global ocean model used in this study is described in detail in Dijkstra and Weijer (2005) to which the reader is referred for full details. In the AMOC model hierarchy (Dijkstra, 2024), this model is located between idealised multi-basin ocean-only models and EMICs. The model captures global ocean flows at low resolution but has a strongly simplified atmospheric model. Through the implicit setup, steady states can be determined versus parameters without using any time stepping, as explained in more detail in the next subsection below.

The governing equations of the ocean model are the hydrostatic, primitive equations in spherical coordinates on a global domain which includes continental geometry as well as bottom topography. The ocean velocities in eastward (zonal) and northward (meridional) directions are indicated by u and v, the vertical velocity is indicated by w, the pressure by p and the temperature and salinity by T and S, respectively. The horizontal resolution of the model is about 4° (a 96 × 38 Arakawa C-grid on a domain [180° W,180°E] × [85.5°S, 85.5°N]) and the grid has 12 vertical levels. The vertical grid is non-equidistant with the surface (bottom) layer having a thickness of 50 m (1000 m), respectively.

Vertical and horizontal mixing of momentum and of tracers (i.e., heat and salt) are represented by a Laplacian formulation with prescribed ‘eddy’ viscosities A_H and A_V and eddy diffusivities K_H and K_V, respectively. As in Dijkstra (2007), we will use the depth dependent values of K_V and K_H (Bryan and Lewis, 1979; England, 1993) given by

(1a)

K V (z) = K V 0 − A s arctan (λ V (z − z *)),

M1 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {K_V}(z) = K_V^0 - {A_s}\arctan ({\lambda _V}(z - {z_*})), \] \end{document}

(1b)

K H (z) = K H 0 + (A r − K H 0) e z λ H,

M2 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {K_H}(z) = K_H^0 + ({A_r} - K_H^0){e^{\frac{z}{{{\lambda _H}}}}}, \] \end{document}

with z ∈ [–5000, 0] m. Here, $K H 0 = 0.5 × 103$ M3 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ K_H^0 = 0.5 \times {10^3} \] \end{document} m²s^–1, A_r = 1.0 × 10³ m²s^–1, $K V 0 = 8.0 × 10 − 5$ M4 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ K_V^0 = 8.0 \times {10^{-5}} \] \end{document} m²s^–1, A_s = 3.3 × 10^–5 m²s^–1, λ_V = 4.5 × 10^–3 m^–1, λ_H = 5 × 10² m and z_* = – 2.5 × 10³ m. A plot of the vertical structure of K_V and K_H can be found in Figure 1 of Dijkstra (2007). In this way, the vertical diffusivity K_V increases from 0.31 × 10^–4 m²s^–1 at the surface to 1.3 × 10^–4 m²s^–1 near the bottom of the flow domain. The horizontal diffusivity K_H increases monotonically from 0.5 × 10³ m²s^–1 at the bottom of the ocean to 1.0 × 10³ m²s^–1 near the surface.

The ocean flow is forced by the observed annual-mean wind stress as given in Trenberth et al. (1989). The upper ocean is coupled to a simple energy-balance atmospheric model (see Appendix in Dijkstra and Weijer (2005)) in which only the heat transport is modelled. The freshwater flux will be prescribed in each of the results in section 3 and the model has no sea-ice component. Both the neglected moisture transport and sea-ice ocean interactions may affect the results below, but these effects are outside the scope of this paper. The surface forcing is represented as a body forcing over the upper layer. On the continental boundaries, no-slip conditions are prescribed and the heat- and salt fluxes are zero. At the bottom of the ocean, both the heat and salt fluxes vanish and slip conditions are assumed.

2.2 Methods

The discretised steady equations can be written as a nonlinear algebraic system of equations of the form

(2)

G (x, μ) = 0,

M5 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\bf{G}}({\bf{x}},\mu) = 0, \] \end{document}

where x is the state vector and μ is one of the parameters of the model. For the global ocean model (with a 4° horizontal resolution and 12 layers in the vertical) the dimension of the state space (and of x) is 96 × 38 × 13 × 6 = 284,544; where the number 13 comes from the 12 ocean levels plus the atmospheric energy balance model, and the 6 from the number of unknowns (u,v,w,p,T,S) per grid point.

We use the pseudo-arclength continuation method (Keller, 1977), where the branch of steady solutions versus μ is parametrised by an arclength s. To close the set of equations (because of the new variable s) the arclength is normalised leading to the equations

(3a)

G (x (s), μ (s)) = 0,

M6 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\bf{G}}({\bf{x}}(s),\mu (s)) = 0, \] \end{document}

(3b)

x ˙ 0 T (x (s) − x 0) + μ ˙ 0 (μ (s) − μ 0) − (s − s 0) = 0,

M7 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\bf{\dot x}}_0^T({\bf{x}}(s) - {{\bf{x}}_0}) + {\dot \mu _0}(\mu (s) - {\mu _0}) - (s - {s_0}) = 0, \] \end{document}

where (x₀, μ₀) is a previously computed solution and the dot indicates differentiation to s. The equation (3b) represents a normalisation of the tangent on the branch of steady solutions. The linear stability of each steady state is determined by solving a generalised eigenvalue problem using the Jacobi-Davidson QZ method (Dijkstra, 2005).

The procedure to compute bifurcation diagrams of the model, including biases in the freshwater forcing, is the following:

Under restoring conditions for the surface salinity field (Levitus, 1994), a steady solution is determined for standard values of the parameters of the model (Dijkstra and Weijer, 2005). From this steady solution the freshwater flux, below referred to the Levitus flux $F S L$ M8 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ F_S^L \] \end{document} , is diagnosed.
A freshwater flux over a region near Newfoundland (Figure 2) with domain [60°W, 24°W] × [54°N, 66°N] is prescribed (in addition to $F S L$ M9 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ F_S^L \] \end{document} ) with strength $γ A F S A$ M10 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\gamma _A}F_S^A \] \end{document} Sv, where $F S A = 1$ M11 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ F_S^A = 1 \] \end{document} in this domain and zero outside. Similarly, a bias freshwater flux is prescribed over the Indian Ocean domain [52°E, 104°E] × [20°S, 10°N] (same red-outlined region as in Figure 1) with amplitude $γ I F S I$ M12 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\gamma _I}F_S^I \] \end{document} Sv, where $F S I = 1$ M13 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ F_S^I = 1 \] \end{document} in this domain and zero outside. The total freshwater flux is prescribed as

(4)

F S = F S L + γ A F S A + γ I F S I − Q F S C,

M14 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {F_S} = F_S^L + {\gamma _A}F_S^A + {\gamma _I}F_S^I - QF_S^C, \] \end{document}

where $F S C = 1$ M15 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ F_S^C = 1 \] \end{document} in a compensation domain (specified below) and the quantity Q is determined such that

(5)

∫ S o a F S r 02 cos θ d θ d ϕ = 0,

M16 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \int_{{S_{oa}}} {{F_S}} r_0^2\cos \theta d\theta d\phi = 0, \] \end{document}

where S_oa is the total ocean surface and r₀ the radius of the Earth.

(iii) In the results below, we will consider two cases of compensation: (a) global compensation, i.e. C is the global ocean domain (Figure 2a) as in (Dijkstra 2007) and (b) C is the Pacific domain (Figure 2b), so there is no compensation over the Atlantic. In each case, for different (but fixed) values of γ_I, a branch of steady solutions versus γ_A is calculated under the freshwater forcing (4), starting from the solutions determined under (i) for γ_A = γ_I = 0.

Figure 2

Areas where freshwater flux anomalies are applied with their strengths γ_A and γ_I; also the global compensation region (a) and the Pacific region (b) is shown.

The Pacific compensation is here purely introduced because there needs to be a compensation surface flux to keep the salinity constant. It is not introduced here to represent the origin of the bias in the Indian Ocean surface freshwater flux or to represent that the excess evaporation of the Pacific water will lead to excess precipitation over the Indian Ocean. Furthermore, γ_A is purely considered as a control parameter to freshen the North Atlantic and not to represent the surface freshwater bias in this region. The global compensation case is discussed here in detail because it has been frequently used in quasi-equilibrium studies (Rahmstorf et al., 2005), including the recent study of van Westen et al. (2024).

3 Results

In the results below, we concentrate on the bifurcation diagrams and freshwater and salt balances. Plots of the typical AMOC patterns for slightly different parameter values can be found in Dijkstra (2007) and no solutions with Pacific sinking are found.

3.1 Global Compensation

The bifurcation diagram for the global compensation case (Figure 2a) with γ_I = 0 (no Indian Ocean freshwater flux bias), where the maximum AMOC strength below 1000 m (Ψ_A) is plotted versus γ_A (both in Sv), is shown as the black curve in Figure 3a. With increasing γ_A, stable (black solid curves) steady states exist for which the AMOC strength decreases and at $γ A 1 = 0.186$ M17 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \gamma _A^1 = 0.186 \] \end{document} Sv, a first saddle-node bifurcation L₁ occurs. With decreasing γ_A, a branch of unstable steady states (black dashed curves) exists down to a second saddle-node bifurcation L₂ at $γ A 2 = 0.054$ M18 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \gamma _A^2 = 0.054 \] \end{document} Sv.

Figure 3

Bifurcation diagrams for the case of global compensation where the maximum strength of the AMOC below 1000 m (Ψ_A) is plotted versus the strength of the anomalous freshwater forcing γ_A for different values of γ_I. Solid (dashed) curves indicate stable (unstable) branches. The dots indicate the saddle-node bifurcations.

The width of the multi-stable regime, often called the hysteresis width Δ_H, is given by

(6)

Δ H = | γ A 1 − γ A 2 | .

M19 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\Delta _H} = |\gamma _A^1 - \gamma _A^2|. \] \end{document}

For the case γ_I = 0, we find Δ_H = 0.132 Sv in this model. In typical quasi-equilibrium model studies (Rahmstorf et al., 2005), where γ_A is varied with about 0.05 Sv/1000 years, the width is typically overestimated. With continuation methods, as used here, one is able to determine the hysteresis width very accurately as the values of $γ A 1, 2$ M20 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \gamma _A^{1,2} \] \end{document} are computed explicitly.

With increasing values of γ_I (adding fresh water over the Indian Ocean) both saddle node-bifurcations L₁ and L₂ move to larger values of γ_A (Figure 3) indicating that the multi-stable regime occurs for higher anomalous additional Atlantic freshwater forcing. Hence, GCMs having a ‘Levitus-like’ surface salinity with the given Indian Ocean freshwater flux bias are less likely to be in a multi-stable regime than models without such a bias. The width of the multi-stable regime versus γ_I does not change much for these γ_I values; Δ_H = 0.137 Sv and Δ_H = 0.138 Sv are found for γ_I = 0.37 Sv and γ_I = 0.52 Sv, respectively. These imposed γ_I values are larger than the CESM biases of about 0.05 Sv (Figure 1). Note that freshwater biases outside the Indian Ocean freshwater region also contribute to a fresher Indian Ocean, in particular the positive biases over land which enhance river run-off.

To understand the shift of the branches in the bifurcation diagram, we consider the Atlantic freshwater transport by the AMOC and the gyres. Following De Vries and Weber (2005), the quantities F_ov (the overturning component) and F_az (the azonal component) are computed as

(7a)

F ov (θ) = − r 0 S 0 ∫ S θ v ¯ (< S > − S 0) d z;

M21 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {F_{{\rm{ov}}}}(\theta) = - \frac{{{r_0}}}{{{S_0}}}\int_{{S_\theta }} {\bar v} (< S > - {S_0})dz; \] \end{document}

(7b)

F az (θ) = − r 0 S 0 ∫ S θ v ′ S ′ ¯ d z .

M22 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {F_{{\rm{az}}}}(\theta) = - \frac{{{r_0}}}{{{S_0}}}\int_{{S_\theta }} {\overline {v^{\prime}S^{\prime}}} dz. \] \end{document}

where S₀ = 35 psu is a reference salinity, r₀ is the radius of the Earth, and S_θ is the boundary (longitude, depth) at latitude θ. Here, the quantities $v ¯, S ¯, < v >$ M23 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \bar v,\bar S, < v > \] \end{document} and <S> are given by

(8a)

v ¯ = ∫ v cos θ d ϕ; < v > = v ¯ ∫ cos θ d ϕ,

M24 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \bar v = \int v \cos \theta d\phi ; < v > = \frac{{\bar v}}{{\int {\cos \,} \theta d\phi }}, \] \end{document}

(8b)

S ¯ = ∫ S cos θ d ϕ; < S > = S ¯ ∫ cos θ d ϕ,

M25 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \bar S = \int S \cos \theta d\phi ; < S > = \frac{{\bar S}}{{\int {\cos } \,\theta d\phi }}, \] \end{document}

and v′ = v – <v> and S′ = S – <S>. The physical meaning of these quantities is extensively discussed in De Vries and Weber (2005) and Dijkstra (2007).

The existence of the saddle-node bifurcations L₁ and L₂ can be connected to the behavior of F_ovS = F_ov (35°S). In simple box models (Cessi, 1994; Rahmstorf, 1996), the saddle-node bifurcation L₁ at $γ A 1$ M26 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \gamma _A^1 \] \end{document} is related to a minimum in F_ovS. In our model, this is more complicated as there is also a gyre-driven freshwater transport, and the F_ovS minimum is only approximate. The saddle-node bifurcation L₂ at $γ A 2$ M27 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \gamma _A^2 \] \end{document} is near to a zero of F_ovS along the upper branch of the AMOC, as discussed at length in Dijkstra (2007). So to explain the shift in positions in the saddle-node bifurcations we focus on the behaviour of F_ovS, F_ovN = F_ov(60°N), while also monitoring their difference ΔF_ov = F_ovS – F_ovN and the associated behaviour of the freshwater transports by the gyres.

The results for F_ovS, F_ovN and ΔF_ov are shown in Figure 4a, with the case γ_I = 0 (γ_I = 0.37 Sv) as solid (dashed) curves. For the chosen northern latitude (60°N), the freshwater flux F_ovN is negative and the AMOC transports freshwater southwards for all values of γ_A. For γ_A = 0, the AMOC transports freshwater northwards at 35°S as F_ovS > 0. With increasing γ_A, F_ovS decreases and, for γ_I = 0, becomes negative close to the saddle-node bifurcation L₂ at $γ A 2$ M28 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \gamma _A^2 \] \end{document} (indicated by the thin vertical line). In other words, the (negative) F_ovS sign indicates the multi-stable regime.

Figure 4

Global compensation case. (a) Values of the AMOC induced freshwater transport at the southern boundary (35°S) of the Atlantic (F_ovS), the northern boundary (60°N) of the Atlantic (F_ovN), and their difference ΔF_ov. (b) Same as (a) but for the azonal transport (F_azS, F_azN and ΔF_az). The vertical thin lines indicate the saddle-node bifurcation L₂. The solid curves are for γ_I = 0, corresponding to the black curve in Figure 3. The dashed curves are for γ_I = 0.37 Sv, corresponding to the blue curve in Figure 3.

For γ_I = 0.37 Sv, the value of F_ovS at γ_A = 0 is much larger than that for γ_I = 0. Because F_ovS decreases with γ_A at the same rate as for γ_I = 0, the location where F_ovS ≈ 0 and hence the bifurcation L₂ occurs at a larger value of γ_A. Also the location where F_ovS obtains its minimum, and hence the position of L₁ (at $γ A 1$ M29 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \gamma _A^1 \] \end{document} ) shifts to the right. The gyre transport changes F_azS, F_azN and ΔF_az with γ_A are shown in Figure 4b and do not change much with γ_I.

As shown in Dijkstra (2007), the fully-implicit model allows for a closed salt balance over the Atlantic from θ_s =35°S to θ_n = 60°N from which changes in advective, diffusive and surface contributions can be determined. The terms in this balance are shown in Figure 5 for both cases γ_I = 0 and γ_I = 0.37 Sv. Expressions for these terms were presented in Dijkstra (2007), but are repeated here for convenience, i.e.

(9a)

Φ s = ∫ S o a S 0 F S r 02 cos θ d ϕ d θ,

M30 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\Phi ^s} = \int_{{S_{oa}}} {{S_0}} {F_S}r_0^2\cos \,\theta d\phi d\theta, \] \end{document}

(9b)

Φ a (θ) = − ∫ S θ v S r 0 cos θ d ϕ d z,

M31 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\Phi ^a}(\theta) = - \int_{{S_\theta }} v S{r_0}\cos \theta d\phi dz, \] \end{document}

(9c)

Φ d (θ) = ∫ S θ K H ∂ S ∂ θ cos θ d ϕ d z,

M32 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\Phi ^d}(\theta) = \int_{{S_\theta }} {{K_H}} \frac{{\partial S}}{{\partial \theta }}\cos \theta d\phi dz, \] \end{document}

where $Φ a$ M33 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\Phi ^a} \] \end{document} and $Φ d$ M34 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\Phi ^d} \] \end{document} are the advective and diffusive fluxes through the boundary S_θ respectively. The overall balance is given by

(10)

Φ b = Φ a (θ n) − Φ d (θ n) − Φ a (θ s) + Φ d (θ s) − Φ s

M35 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\Phi ^b} = {\Phi ^a}({\theta _n}) - {\Phi ^d}({\theta _n}) - {\Phi ^a}({\theta _s}) + {\Phi ^d}({\theta _s}) - {\Phi ^s} \] \end{document}

Figure 5

Global compensation case. Terms in the integrated salt balance over the Atlantic basin over the upper branches in Figure 3 (up to $γ A 1$ M36 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \gamma _A^1 \] \end{document} ), with expressions of the terms as indicated in equation (9) and equation (10) for θn = 60°N and θs = 35°S. The solid curves are for the case γI = 0 and the dashed ones for γI = 0.37 Sv.

Indeed, the term Φ^b is much smaller than the individual terms (black curves in Figure 5) giving a nearly closed salt balance over the Atlantic basin for all values of the parameters.

First consider the case γ_I = 0 (solid curves in Figure 5) and the upper branch in the bifurcation diagram up to L₁. For γ_A = 0, the surface (virtual) salt flux is approximately balanced by the fluxes at the southern boundary. The surface evaporation is larger than the precipitation (Φ^s > 0) and this salt is transported out of the Atlantic basin at the southern boundary (Φ^a(θ_s) < 0). The fact that the diffusive flux is relatively large here (compared to typical GCMs) is that the model has a coarse resolution so needs a relatively high horizontal diffusivity to prevent wiggles to occur near the boundaries. The salt fluxes at the northern boundary are less important and the respective components are about a factor 2 to 4 smaller than at the southern boundary.

With increasing γ_A, the surface salt flux decreases as freshwater is put into the North Atlantic. The diffusive salt transports do not respond but the southward advective salt transport at θ_s weakens. As the diffusive flux is directed to transport salt into the basin at the southern boundary, the value of γ_A where the sign change in salt transport occurs is around 0.1 Sv (Figure 5). Note that the gyre and AMOC components cannot be distinguished in the advective fluxes.

When γ_I = 0.37 Sv (dashed curves in Figure 5), the surface salt flux increases for γ_A = 0 compared to the case γ_I = 0. This is due to the global compensation as a negative salt flux in the Indian Ocean is compensated by a positive one over part of the Atlantic. Hence, the curve for Φ^s shifts upwards and so the compensating advective flux at the southern boundary shifts downwards. A second effect is that the changed surface freshwater flux pattern leads to a modified salinity distribution in the Atlantic. This increases the AMOC (Figure 3) strength and hence also its salt transport out of the basin at the southern boundary.

Because the diffusive fluxes are not much affected by the Indian Ocean freshwater input, the fluxes Φ^s and Φ^a(θ_s) change with γ_A in the same way as for the case γ_I = 0. Similarly, the starting value of Φ^s at γ_A = 0 is now larger and it takes a larger value of γ_A to change the sign of the freshwater flux at the southern boundary and to reach a minimum in this quantity. Hence the saddle-node bifurcations L₁ and L₂ shift to larger values of γ_A.

3.2 Pacific Compensation

Since the global compensation has a substantial influence on the position of the saddle-node bifurcations (Figure 3), we now consider the case where compensation is only over the Pacific domain as indicated in Figure 2b. The bifurcation diagrams in this case are, for different values of γ_I, shown in Figure 6. The shift of the saddle-nodes to larger values of γ_A is much smaller than for the global compensation case (Figure 3). The hysteresis width itself for γ_I = 0 with Δ_H = 0.103 Sv for the Pacific compensation case is a bit smaller than for the global compensation case (Δ_H = 0.138 Sv). This width is only slightly larger for the case γ_I = 0.52 Sv, i.e. Δ_H = 0.111 Sv.

Figure 6

Bifurcation diagram for the case of Pacific compensation where the maximum strength of the AMOC below 1000 m (Ψ_A) is plotted versus the strength of the anomalous Atlantic freshwater forcing γ_A for different values of γI. The dots indicate the saddle-node bifurcations.

For the analysis of the freshwater and salt balances, we choose the larger value (γ_I = 0.52 Sv) instead of γ_I = 0.37 Sv (used in the global compensation case), as the differences are more clearly visible. The freshwater transports by the AMOC and by the gyres (Figure 7) show that for γ_A = 0, F_ovS is larger for γ_I = 0.52 Sv (Figure 7a), compared to the case γ_I = 0.0 Sv. Hence also changes in the freshwater balance are induced in the Atlantic, but both the direct effect of compensation of the surface freshwater flux and the secondary effect of an AMOC increase are much smaller. Of course, in such a diffusive and viscous model, the Agulhas retroflection is in a diffusive retroflection regime (Dijkstra and De Ruijter, 2001) and there is additional freshwater transport from the Indian to the Atlantic when γ_I > 0. As the wind-driven freshwater transport does not change much with γ_A (Figure 7b), the AMOC transports more freshwater into the basin and hence F_ovS becomes more positive.

Figure 7

Pacific compensation case. (a) Values of the AMOC induced freshwater transport at the southern boundary (35°S) of the Atlantic (F_ovS), the northern boundary (60°N) of the Atlantic (F_ovN), and their difference (ΔF_ov). (b) Same as (a) but for the azonal transport (F_azS, F_azNand ΔF_az). The vertical thin lines indicate the saddle-node bifurcation L₂. The solid curves are for γ_I = 0, corresponding to the black curve in Figure 6. The dashed curves are for γ_I = 0.52 Sv, corresponding to the red curves in Figure 6.

In the Atlantic salt balance (Figure 8) the surface salt flux is indeed the same for both values of γ_I, because there is no compensation anymore over the Atlantic. The advective transport of salt becomes slightly more negative over the southern boundary for γ_I = 0.52 Sv indicating that indeed a small amount of freshwater is transported into the Atlantic basin. Also the diffusive transport of salt (into the basin) decreases and this approximately balances the advective contribution. As the advective contribution is much smaller than the compensation contribution in the global compensation case, the shift of the saddle-node bifurcations to larger values of γ_A is much smaller.

Figure 8

Pacific compensation case. Terms in the integrated salt balance over the Atlantic basin over the upper branch in Figure 6, with expressions of the terms as indicated in the main text. The solid curves are for the case γ_I = 0 and the dashed ones for γ_I = 0.52 Sv.

4 Summary and Discussion

Following earlier studies on CMIP3 and CMIP5 models (Drijfhout et al., 2013; Mecking et al., 2016), also many CMIP6 models have large biases in surface freshwater fluxes which lead to an AMOC with an Atlantic freshwater transport that is in disagreement with observations (van Westen and Dijkstra, 2024). The most important model bias is in the Atlantic Surface Water properties, which arises from the surface freshwater flux over the Indian Ocean giving too fresh water entering the Atlantic through the Agulhas Leakage.

In this paper, we have addressed the effects of this surface freshwater flux bias on the multiple equilibrium regime of the AMOC using the fully implicit global ocean-atmosphere model of Dijkstra and Weijer (2005), for which explicit bifurcation diagrams can be computed. This is a fantastic capability as both stable and unstable steady states can be computed, and the width of the multiple equilibrium regime can be determined accurately. However, this can only be done with quite a simplified global model, with relatively low resolution and hence being more viscous and diffusive than state-of-the-art, even low-resolution, ocean models. The atmosphere model is only an energy balance model with a prescribed freshwater forcing. The effects of lower diffusivity/viscosity can be well anticipated as further instabilities will occur (e.g. baroclinic/barotropic) of the steady states. The effect of an active moisture flux in the atmosphere will lead to changing freshwater forcing with a varying AMOC which can also lead to shifts in the equilibria (Toom et al., 2012). In general, in terms of quantitative results on changes of the bifurcation diagrams, the model is probably not that useful.

Qualitatively, however, the model provides very useful information, as it indicates that the multiple equilibrium regime does not disappear due to the Indian Ocean freshwater flux biases, but that it shifts to higher values of the North Atlantic anomalous freshwater flux. This shift is dependent on the way the latter flux is compensated. Surface salinity patterns in the Atlantic depend on the compensation of the Indian Ocean bias and, in both compensation cases considered here, lead to a slight increase in the AMOC. In the Levitus background state (with F_ovS > 0), this leads to a larger transport of salt out of the Atlantic basin. Hence, a larger anomalous North Atlantic surface freshwater flux is needed to activate the basin-wide salt advection feedback, i.e., to a situation where the AMOC exports fresh water. When there is compensation of the Indian Ocean freshwater flux in the Atlantic, the Atlantic becomes saltier and hence the AMOC transports more salt out of the basin. This leads to an additional, and here larger, shift of the bifurcation diagram compared to when there is no compensation in the Atlantic.

The mechanisms identified are useful for interpreting results from GCMs and for designing new simulations with these models. First it shows that positive biases in surface freshwater flux lead to shifts in the bifurcation diagram which would imply that such multiple equilibrium regimes (and hence AMOC collapses) could exist in these models but are located in a parameter regime where one would not normally perform simulations. Second, it provides a hint why efforts to find AMOC collapses in these models may not have been successful (Mecking et al., 2016; Jackson and Wood, 2018b, a). Even an enormous freshwater input in a parameter regime which is outside the multiple equilibrium regime would only lead to a weakened (but no collapsed) AMOC.

To obtain a better confidence in AMOC stability, the identification of whether the AMOC is in a multiple equilibrium regime is crucial. Our study suggests several ways forward to find an AMOC collapsed state in state-of-the-art models. One may perform a long quasi-equilibrium simulation up to very large freshwater flux input such as in the FAMOUS model (Hawkins et al., 2011) to find the collapse (this may take a few thousand years of simulation, so is expensive). In doing this, it is better to compensate outside of the Atlantic when using surface fluxes, as the latter will introduce an additional shift and so one has to integrate longer to find an AMOC collapse. Such compensation procedures (including compensation over the volume) are now also more common in GCMs (Jackson et al., 2022). The alternative is to address and improve the biases in the atmospheric components of the models, as illustrated here in the case of the Indian Ocean bias, but this is not an easy issue. The origin of these biases may even be a coupled problem, as the bias strength is positively correlated with the AMOC strength (van Westen and Dijkstra, 2024).

We hope that this study will motivate the design of new simulations in state-of-the-art GCMs to assess the multiple equilibrium regime of the AMOC. Detection of such a multiple equilibrium regime would have a large impact on climate change research and probability estimates of AMOC tipping under global warming (Armstrong McKay et al., 2022) would likely need to be revised.

Acknowledgements

The authors thank Dr. Fred Wubs (University of Groningen, NL) for helping with reviving the old Fortran code of the model (Dijkstra, 2007) used in this paper. Model output for the CESM simulations can accessed at https://ihesp.github.io/archive/ and the ERA5 reanalysis product is available at https://doi.org/10.24381/cds.f17050d7.

Funding Information

H.A.D. and R.M.v.W. are funded by the European Research Council through the ERC-AdG project TAOC (project 101055096).

Competing Interests

The authors have no competing interests to declare.

References

Armstrong McKay, DI, Staal, A, Abrams, JF, Winkelmann, R, Sakschewski, B, Loriani, S, Fetzer, I, Cornell, SE, Rockström, J and Lenton, TM. 2022. Exceeding 1.5 C global warming could trigger multiple climate tipping points. Science, 377(6611): eabn7950. DOI: https://doi.org/10.1126/science.abn7950
Bryan, K and Lewis, LJ. 1979. A water mass model of the world ocean. J. Geophys. Res., 84: 2503–2517. DOI: https://doi.org/10.1029/JC084iC05p02503
Bryden, HL, King, BA and McCarthy, GD. 2011. South Atlantic overturning circulation at 24S. Journal of Marine Research, 69: 38–55. DOI: https://doi.org/10.1357/002224011798147633
Cessi, P. 1994. A simple box model of stochastically forced thermohaline flow. J. Phys. Oceanogr., 24: 1911–1920. DOI: https://doi.org/10.1175/1520-0485(1994)024<1911:ASBMOS>2.0.CO;2
Cimatoribus, A, Drijfhout, S, Toom, den M and Dijkstra, HA. 2012. Sensitivity of the Atlantic meridional overturning circulation to South Atlantic freshwater anomalies. Climate Dynamics, 39: 2291–2306. DOI: https://doi.org/10.1007/s00382-012-1292-5
Cini, M, Zappa, G, Ragone, F and Corti, S. 2024. Simulating AMOC tipping driven by internal climate variability with a rare event algorithm. npj Climate and Atmospheric Science, 7(1): 31. DOI: https://doi.org/10.1038/s41612-024-00568-7
De Niet, A, Wubs, F, van Scheltinga, AT and Dijkstra, HA. 2007. A tailored solver for bifurcation analysis of ocean-climate models. Journal of Computational Physics, 227(1): 654–679. DOI: https://doi.org/10.1016/j.jcp.2007.08.006
De Vries, P and Weber, SL. 2005. The Atlantic freshwater budget as a diagnostic for the existence of a stable shut down of the meridional overturning circulation. Geophys. Res. Letters, 32: L09606. DOI: https://doi.org/10.1029/2004GL021450
Dijkstra, HA. 2005. Nonlinear Physical Oceanography: A Dynamical Systems Approach to the Large Scale Ocean Circulation and El Niño, 2nd Revised and Enlarged edition. Springer, New York, 532 pp.
Dijkstra, HA. 2007. Characterization of the multiple equilibria regime in a global ocean model. Tellus, 59A: 695–705. DOI: https://doi.org/10.1111/j.1600-0870.2007.00267.x
Dijkstra, HA. 2024. The role of conceptual models in climate research. Physica D: Nonlinear Phenomena, 457: 133984. DOI: https://doi.org/10.1016/j.physd.2023.133984
Dijkstra, HA and De Ruijter, WPM. 2001. On the physics of the Agulhas current: Steady retroflection regimes. J. Phys. Oceanogr., 31: 2971–2985. DOI: https://doi.org/10.1175/1520-0485(2001)031<2971:OTPOTA>2.0.CO;2
Dijkstra, HA and Weijer, W. 2005. Stability of the global ocean circulation: basic bifurcation diagrams. J. Phys. Oceanogr., 35: 933–948. DOI: https://doi.org/10.1175/JPO2726.1
Drijfhout, S, Gleeson, E, Dijkstra, HA and Levina, V. 2013. Spontaneous abrupt climate change due to an atmospheric blocking– sea-ice–ocean feedback in an unforced climate model simulation. Proceedings National Acad. Sciences. DOI: https://doi.org/10.1073/pnas.1304912110
England, MH. 1993. Representing the global-scale water masses in ocean general circulations models. J. Phys. Oceanogr., 23: 1523–1552. DOI: https://doi.org/10.1175/1520-0485(1993)023<1523:RTGSWM>2.0.CO;2
Gent, PR. 2018. A commentary on the atlantic meridional overturning circulation stability in climate models. Ocean Modelling, 122: 57–66. DOI: https://doi.org/10.1016/j.ocemod.2017.12.006
Hawkins, E, Smith, RS, Allison, LC, Gregory, JM, Woollings, TJ, Pohlmann, H and De Cuevas, B. 2011. Bistability of the Atlantic overturning circulation in a global climate model and links to ocean freshwater transport. Geophysical Research Letters, 38(10): L10605. DOI: https://doi.org/10.1029/2011GL047208
Hu, A, Meehl, GA, Han, W, Timmermann, A, Otto-Bliesner, B, Liu, Z, Washington, WM, Large, W, Abe-Ouchi, A, Kimoto, M, Lambeck, K and Wu, B. 2012. Role of the Bering Strait on the hysteresis of the ocean conveyor belt circulation and glacial climate stability. Proceedings of the National Academy of Sciences, 109(17): 6417–6422. DOI: https://doi.org/10.1073/pnas.1116014109
Huisman, SE, den Toom, M, Dijkstra, HA and Drijfhout, S. 2010. An Indicator of the Multiple Equilibria Regime of the Atlantic Meridional Overturning Circulation. Journal Of Physical Oceanography, 40(3): 551–567. DOI: https://doi.org/10.1175/2009JPO4215.1
Jackson, LC, Asenjo, EAd, Bellomo, K, Danabasoglu, G, Haak, H, Hu, A, Jungclaus, J, Lee, W, Meccia, VL, Saenko, O, Shao, A, and Swingedouw, D. 2022. Understanding AMOC stability: the North Atlantic Hosing Model Intercomparison Project. Geoscientific Model Development Discussions, 1–32. DOI: https://doi.org/10.5194/gmd-2022-277
Jackson, LC, Kahana, R, Graham, T, Ringer, MA, Woollings, T, Mecking, JV and Wood, RA. 2015. Global and European climate impacts of a slowdown of the AMOC in a high resolution GCM. Climate Dynamics, 45(11): 3299–3316. DOI: https://doi.org/10.1007/s00382-015-2540-2
Jackson, LC and Wood, RA. 2018a. Hysteresis and Resilience of the AMOC in an Eddy-Permitting GCM. Geophysical Research Letters, 45(16): 8547–8556. DOI: https://doi.org/10.1029/2018GL078104
Jackson, LC and Wood, RA. 2018b. Timescales of AMOC decline in response to fresh water forcing. Climate Dynamics, 51(4): 1333–1350. DOI: https://doi.org/10.1007/s00382-017-3957-6
Keller, HB. 1977. Numerical solution of bifurcation and nonlinear eigenvalue problems. In Rabinowitz, PH (ed.), Applications of Bifurcation Theory. Academic Press, New York, U.S.A.
Lenton, TM, Held, H, Kriegler, E, Hall, JW, Lucht, W, Rahmstorf, S and Schellnhuber, HJ. 2008. Tipping elements in the Earth’s climate system. Proceedings of the National Academy of Sciences of the United States of America, 105(6): 1786–93. DOI: https://doi.org/10.1073/pnas.0705414105
Levitus, S. 1994. World Ocean Atlas 1994, Volume 4: Temperature. NOAA/NESDIS E, US Department of Commerce, Washington DC, OC21, 1–117.
Liu, W, Fedorov, AV, Xie, S-P and Hu, S. 2020. Climate impacts of a weakened Atlantic Meridional Overturning Circulation in a warming climate. Science Advances, 6(26): eaaz4876. DOI: https://doi.org/10.1126/sciadv.aaz4876
Liu, W, Xie, S-P, Liu, Z and Zhu, J. 2017. Overlooked possibility of a collapsed Atlantic Meridional Overturning Circulation in warming climate. Science Advances, 3(1): e1601666. DOI: https://doi.org/10.1126/sciadv.1601666
Lohmann, J, Castellana, D, Ditlevsen, PD and Dijkstra, HA. 2021. Abrupt climate change as rate-dependent cascading tipping point. Earth System Dynamics Discussions, 2021, 1–25. DOI: https://doi.org/10.5194/esd-2021-7
Lohmann, J, Dijkstra, HA, Jochum, M, Lucarini, V and Ditlevsen, PD. 2024. Multistability and intermediate tipping of the Atlantic Ocean circulation. Science Advances, 10(12): eadi4253. DOI: https://doi.org/10.1126/sciadv.adi4253
Marotzke, J. 2000. Abrupt climate change and thermohaline circulation: Mechanisms and predictability. Proc. Natl. Acad. Sci., 97: 1347–1350. DOI: https://doi.org/10.1073/pnas.97.4.1347
Mecking, JV, Drijfhout, SS, Jackson, LC and Graham, T. 2016. Stable AMOC off state in an eddy-permitting coupled climate model. Climate Dynamics, 47(7): 2455–2470. DOI: https://doi.org/10.1007/s00382-016-2975-0
Mulder, TE, Goelzer, H, Wubs, FW and Dijkstra, HA. 2021. Snowball Earth Bifurcations in a Fully-Implicit Earth System Model. International Journal of Bifurcation and Chaos, 31(06): 2130017. DOI: https://doi.org/10.1142/S0218127421300172
Orihuela-Pinto, B, England, MH and Taschetto, AS. 2022. Interbasin and interhemispheric impacts of a collapsed atlantic overturning circulation. Nature Climate Change, 12(6): 558–565. DOI: https://doi.org/10.1038/s41558-022-01380-y
Rahmstorf, S. 1996. On the freshwater forcing and transport of the Atlantic thermohaline circulation. Climate Dynamics, 12(12): 799–811. DOI: https://doi.org/10.1007/s003820050144
Rahmstorf, S, Crucifix, M, Ganopolski, A, Goosse, H, Kamenkovich, I, Knutti, R, Lohmann, G, March, R, Mysak, L, Wang, Z and Weaver, AJ. 2005. Thermohaline circulation hysteresis: a model intercomparison. Geophys. Res. Letters, L23605.
Stommel, H. 1961. Thermohaline convection with two stable regimes of flow. Tellus, 2: 244–230. DOI: https://doi.org/10.1111/j.2153-3490.1961.tb00079.x
Stouffer, RJ, Yin, J, Gregory, JM, Dixon, KW, Spelman, MJ, Hurlin, W, Weaver, AJ, Eby, M, Flato, GM, Hasumi, H, Hu, A, Jungclaus, JH, Kamenkovich, IV, Levermann, A, Montoya, M, Murakami, S, Nawrath, S, Oka, A, Peltier, WR, Robitaille, DY, Sokolov, AP, Vettoretti, G and Weber, SL. 2006. Investigating the causes of the response of the thermohaline circulation to past and future climate changes. Journal of Climate, 19: 1365–1387. DOI: https://doi.org/10.1175/JCLI3689.1
Toom, den M, Dijkstra, HA, Cimatoribus, AA and Drijfhout, SS. 2012. Effect of Atmospheric Feedbacks on the Stability of the Atlantic Meridional Overturning Circulation. Journal of Climate, 25(12): 4081–4096. DOI: https://doi.org/10.1175/JCLI-D-11-00467.1
Trenberth, KE, Olson, JG and Large, WG. 1989. A global ocean wind stress climatology based on ECMWF analyses. Technical report, National Center for Atmospheric Research, Boulder, CO, U.S.A.
van Westen, RM and Dijkstra, HA. 2023. Asymmetry of AMOC Hysteresis in a State-Of-The-Art Global Climate Model. Geophysical Research Letters, 50(22): e2023GL106088. DOI: https://doi.org/10.1029/2023GL106088
van Westen, RM and Dijkstra, HA. 2024. Persistent Climate Model Biases in the Atlantic Ocean’s Freshwater Transport. Ocean Science, 20: 549–567. DOI: https://doi.org/10.5194/os-20-549-2024
van Westen, RM, Kliphuis, M and Dijkstra, HA. 2024. Physics-based early warning signal shows that AMOC is on tipping course. Science Advances, 10(6): eadk1189. DOI: https://doi.org/10.1126/sciadv.adk1189
Vellinga, M, Wood, RA and Gregory, JM. 2002. Processes governing the recovery of a perturbed thermohaline circulation in HadCM3. J. Climate, 15: 764–780. DOI: https://doi.org/10.1175/1520-0442(2002)015<0764:PGTROA>2.0.CO;2
Weijer, W, Cheng, W, Drijfhout, SS, Fedorov, AV, Hu, A, Jackson, LC, Liu, W, McDonagh, EL, Mecking, JV and Zhang, J. 2019. Stability of the Atlantic Meridional Overturning Circulation: A Review and Synthesis. Journal Of Geophysical Research-Oceans, 124(8): 5336–5375. DOI: https://doi.org/10.1029/2019JC015083

Tellus A: Dynamic Meteorology and Oceanography