Rui Wang
School of Mathematical Sciences, Fudan University, Shanghai, China.
DOI: 10.4236/ojmsi.2025.134013   PDF    HTML   XML   55 Downloads   341 Views   Citations

Abstract

We use a three-box Stommel model not only to study the expected trajectory of the Atlantic Meridional Overturning Circulation (AMOC), but also to determine the mechanisms underlying Dansgaard-Oeschger (D-O) cycles. By constructing this enhanced model framework, we demonstrate that the inclusion of a third-box accounting for global freshwater influences potentially accelerates the AMOC collapse timeline compared to conventional two-box model predictions. Regarding D-O cycles, our analysis identifies multidecadal variability in AMOC dynamics and the polar see-saw hypothesis (where hemispheric glaciation/destruction requires prolonged temporal lag to manifest reciprocal effects) as primary mechanisms. Through temporal lag implementation in our three-box model, we successfully reproduced characteristic D-O cycle patterns, proposing a novel explanatory framework for these paleoclimatic oscillations. This modeling approach bridges freshwater forcing, AMOC instability, and abrupt climate transitions during glacial periods.

Keywords

Stommel Model, AMOC, Dansgaard-Oeschger Cycles, Temporal Lag

Share and Cite:

1. Introduction

In this paper, we primarily address two key research questions: the projected trajectory of AMOC and the construction and mechanistic analysis of D-O event models.

1.1. Introduction of AMOC and Dansgaard-Oeschger Cycles

The Atlantic Meridional Overturning Circulation (AMOC) is an ocean current system transporting warm water in the upper 1 km of the Atlantic Ocean from the tropics towards the North Atlantic, where it cools, becomes denser, sinks, and returns southwards at depth between 2 - 3 km at the magnitude of 20 Sverdrups (Sv) [1] (see Figure 1 for sea surface temperatures of the North Atlantic). This circulation regulates global climate by distributing heat, influencing weather patterns and marine ecosystems across the world [2]-[5].

Figure 1. Sea surface temperatures of the North Atlantic [6].

Dansgaard-Oeschger (D-O) events are millennial-scale climatic oscillations (1470 - 3000 year intervals) observed in Greenland ice cores [7]-[10]. These events consist of abrupt spike-like events with temperatures jumping up to 8˚C - 16˚C over several decades, which continues for centuries, followed by slow cooling [11]. Key lines of evidence include isotopic δ18O shifts in ice cores from Greenland, such as GISP2 and matched temporal changes in methane levels across the entire hemisphere, indicating atmospheric responses to large-scale events (115 - 12 ka) (see Figure 2).

Figure 2. Lsotope data for Antarctic and Greenland ice cores.

The polar see-saw is the phenomenon that temperature changes in the northern and southern hemispheres may be out of phase [12]-[14]. The polar see-saw hypothesis posits that when glaciation intensifies or depletes, it requires substantial time to exert impacts on the opposite hemisphere. We consider the polar see-saw hypothesis as one of the crucial factors in the formation of Dansgaard-Oeschger cycles.

1.2. Related Works

Recent research has uncovered swift reductions in the strength of the Atlantic Meridional Overturning Circulation (AMOC) based on both current observations and analyses [15]-[18], and multiple studies utilizing models have assessed the AMOC and predicted that it could collapse [19]-[25]. Many researchers have explored the AMOC using Stommel-type models. In [26], the authors employed the two-box Stommel Model by dividing the North Atlantic into two boxes (high-latitude and low-latitude) and examined certain properties. However, this model neglects the presence of other global oceans, particularly the Pacific. Observational evidence of global thermohaline circulation pathways also reveals interactions between the South Atlantic and the Pacific. Subsequently, Stommel’s three-box model [27], Stommel’s four-box model [28] and Stommel’s six-box model [29] were developed to simulate a self-sustained oscillatory circulation of AMOC.

Many of the climate signals observed across global systems are predominantly attributable to AMOC variability modulated by freshwater forcing anomalies [30]. In [31] [32], the author suggested that the D-O cycles were caused by multiple equilibrium or variations in the strength of AMOC. Subsequently, in [33] [34], the authors suggested that the AMOC changes are promoted by variations in the Arctic ice sheet runoff, while freshwater runs into the North Atlantic might occur in association with changes of the ice margin, respectively. Then, in [35], the authors provide a simple explanation for D-O cycles by employing the freshwater starvation mechanism: a strong AMOC elevates temperatures, increasing freshwater runoff into the Arctic. This progressively freshens the water until the salinity entering the North Atlantic drops low enough to trigger a transition of the AMOC into a weaker state. In this weakened state, colder temperatures over the ice sheet reduce freshwater input to the Arctic, allowing salinity to rise. This eventually drives the AMOC back to a stronger state.

However, there exists an issue with the model used in [35], i.e., it does not consider the large time lag associated with D-O events, such that decades after the sudden AMOC strengthens and causes more freshwater input into the North Atlantic, it takes about a century longer before this larger freshwater input influences the AMOC again. The original publication neglected these time lags. We analyze this aspect in detail in Section 4.

Notably, authors in [36] have introduced temporal lag to establish a self-sustained oscillatory AMOC model. Starting from the deterministic two-box Stommel model, the authors introduced a key modification: the assumption of an approximately 25-year delay in the freshwater input response to AMOC strength changes. This inclusion of a time delay endowed the previously deterministic system with intrinsic oscillatory behavior, mimicking the periodicity of D-O events. Motivated by this, we postulate that analogous mechanisms governing persistent periodicity may exist in D-O cycles.

1.3. Our Contributions

Our paper focuses on two issues: the model simulation of the present-day AMOC and the model simulation of the D-O cycles. Both frameworks are developed upon the three-box Stommel model architecture. Our contributions are listed as follows:

1) In the examination of today’s AMOC, a new box to represent the whole global ocean was added in our exploration, and we took the 3-box Stommel model with the earth-sized third box to probe the AMOC variation. Meanwhile, we verify the stability of the proposed system and apply Matlab to show that the third box makes sufficient influence on the steady state solution of AMOC.

2) In our investigation of D-O cycles, we begin our work by employing the three-box Stommel model to study these cycles. We consider the polar see-saw hypothesis equally vital in the formation of Dansgaard-Oeschger cycles. Using the three-box Stommel model to simulate AMOC, we incorporated temporal lag to replicate this delayed effect, with a characteristic temporal lag of 300 years. Notably, existing models of D-O cycles have universally omitted temporal lag considerations. We successfully reproduced Dansgaard-Oeschger cycles, thereby proposing a new explanatory framework for this phenomenon.

1.4. Outline of Our Paper

The rest of the paper is organized as follows. We will introduce the two-box Stommel Model in Section 2. In Section 3, we will construct a three-box Stommel model pertaining to AMOC and investigate its steady states and changing trends. In Section 4, we will develop a model for Dansgaard-Oeschger (D-O) events based on the three-box Stommel framework to elucidate the underlying mechanisms of D-O events. Sections 3 and 4 constitute parallel components of the study design.

1.5. Parameters and Variables

Here, we introduced some parameters and variables for the present work, which will be employed throughout this paper.

Variable Name	Variable Unit	Variable Meaning
q	Sv (106 × m3/s)	strength of ocean current
$\rho$	kg/m3	density
v	m3	volume of the ocean
T	K	temperature of the ocean
c	K	temperature of the environment
t	year	time
S	ppt	salinity
$F_s$	m/year	freshwater forcing
${\alpha }_t\left(=0.2\right)$	(kg/m3)K−1	constant
${\beta }_s\left(=0.8\right)$	(kg/m3)ppt−1	constant
k	(m3/s)/(kg/m3)	constant
R	m3/s	runoff rate
r	m3/s	speed of heat conducting into water
d	year	temporal lag

Remark: Given the extensive set of variables, we suppress the subscripts of variables in the table. Subscripts will be introduced to variable names in the following text. This notation is adopted for clarity and does not alter the definition of the variables. As an illustration, Section 4 uses $v_1,v_2,v_3$ to represent the volume of the Arctic Ocean, the subpolar North Atlantic, and the subtropical North Atlantic, respectively.

2. Preliminaries

The Stommel two-box model [37] is one of the first mathematical formalism employed to study the dynamical characteristics of AMOC to nonlinear climate forcings. This section references the work in [26] and sketches its main derivations and analytical results relative to the tipping point behaviour.

2.1. The Formulation of Two-Box Stommel Sodel

We divide the Atlantic into two vertically mixed boxes: the high latitudes part of the Atlantic Ocean, the red box, and the low latitudes one. They are connected by AMOC-driven volume transport $q$ . The circulation is governed by the density difference:

$q=k\left({\rho }_1-{\rho }_2\right)=k\left(-{\alpha }_t\text{Δ}T+{\beta }_s\text{Δ}S\right),$ (1)

where $\text{Δ}T=T_1-T_2<0$ (fixed temperature gradient), $\text{Δ}S=S_1-S_2$ (salinity difference), and $k$ is a proportionality constant.

2.2. Salt Budget Equations

Assuming equal box volumes $v_1=v_2=v$ , the salt conservation equations under freshwater forcing $F_s$ become:

$v\frac{\text{d}{S}_1}{\text{d}t}=\left|q\right|\left(S_2-S_1\right)-F_s,$ (2)

$v\frac{\text{d}{S}_2}{\text{d}t}=\left|q\right|\left(S_1-S_2\right)+F_s.$ (3)

Subtracting these yields the governing equation for salinity difference:

$-v\frac{\text{dΔ}S}{\text{d}t}=2\left|q\right|\text{Δ}S+2F_s.$ (4)

Substituting $q$ from Equation (4) and defining rescaled variables $X={\alpha }_T\text{Δ}T$ and $Y={\beta }_R\text{Δ}S$ , we have:

$\left|Y-X\right|Y=-\frac{{\beta }_{s}F}{k}.$ (5)

For $Y<X$ , the equation becomes:

$Y^2-XY-\frac{{\beta }_s{F}_s}{k}=0,$

with one valid solution:

$Y=\frac{X}{2}-\frac{1}{2}{\left(X^2+4\frac{{\beta }_s{F}_s}{k}\right)}^{1/2}.$

For $Y>X$ , the equation becomes:

$Y^2-XY+\frac{{\beta }_s{F}_s}{k}=0,$

with solutions:

$Y=\frac{X}{2}\pm \frac{1}{2}{\left(X^2-4\frac{{\beta }_s{F}_s}{k}\right)}^{1/2}.$

Note that the plus solution for Y in the first case is positive and is therefore not consistent with the assumption $Y<X$ used to obtain that solution, hence the system exhibits up to three equilibria depending on $F_s$ .

2.3. Multiple Equilibria and Stability

Equation (5) has up to three solutions (see Figure 3):

Figure 3. Steady states of salinity difference and AMOC as a function of $F_s$ .

• Strong AMOC (green branch): High $q$ state maintained by salt-advection feedback.

• Unstable AMOC (blue branch): Intermediate state prone to disturbances.

• Collapsed AMOC (red branch): Weak reversed flow under strong freshwater forcing.

Analysis shows that the blue branch solution is unstable: perturbation grows exponentially in this case, whereas deviations from a stable state decrease. It leads to jumps when crossing critical forcing thresholds because of the structure of the bifurcation.

2.4. Tipping Points and Hysteresis

Gradual freshwater forcing increase (e.g., from ice melt under warming) weakens AMOC along the green branch until collapse occurs at $F\approx 1.2$ m/s. Importantly, reducing F to less than the threshold value of about $F\approx 1.2$ m/s is not sufficient to recover the initial state; a lower $F_s$ needs to be sustained in order to obtain this state. The non-reversible nature of this transition results from the quadratic nonlinearity found in Equation (5), typical of abrupt tipping mechanisms in the climate system.

3. A Simple Three-Box Model for AMOC

3.1. Construction of the Three-Box Model

We divide the world ocean into three boxes: the blue box represents the high latitudes part of the Atlantic Ocean, the red box, the low latitudes one and the black box, the global ocean except the Atlantic Ocean. We let the ocean box temperature $T_1,T_2,T_3$ to be constants. Assuming the volumes of the ocean boxes are constant, let $v_1,v_2,v_3$ represent the volume of the blue box, the red box and the black box.

Denote the circulation between the red box and the blue box as $q_1$ , the circulation between the red box and the black box as $q_2$ . The unknowns in the three-box model is the three boxes and the circulation $q_1,q_2$ . See Figure 4.

Figure 4. The three-box Stommel model.

Similar to that in the two-box model, denote the freshing effect caused by glacier melting and excess precipitation in the blue box and the black box as $-F_{s1}$ and $-F_{s3}$ . Considering both the salt transfer between the red box and the blue box, the red box and the black box, we derived

$v_1\frac{\text{d}{S}_1\left(t\right)}{\text{d}t}=\left|q_1\right|\left(S_2-S_1\right)-F_{s1}$

$v_2\frac{\text{d}{S}_2\left(t\right)}{\text{d}t}=\left|q_1\right|\left(S_1-S_2\right)+\left|q_2\right|\left(S_3-S_2\right)+F_{s2}$

$v_3\frac{\text{d}{S}_3\left(t\right)}{\text{d}t}=\left|q_2\right|\left(S_2-S_3\right)-F_{s3}$

Let

$\text{Δ}{T}_1=T_1-T_2<0,\text{Δ}{T}_2=T_3-T_2<0$

$\text{Δ}{S}_1=S_1-S_2,\text{Δ}{S}_2=S_3-S_2$

The circulation is assumed to be proportional to the density difference between the boxes, so $q_i,i=1,2$ can be written as $k_i\left(-{\alpha }_t\text{Δ}{T}_i+{\beta }_s\text{Δ}{S}_i\right),i=1,2$ . Therefore, we derived

$v_1\frac{\text{d}{S}_1\left(t\right)}{\text{d}t}=\left|k_1\left(-{\alpha }_t\left(T_1-T_2\right)+{\beta }_s\left(S_1-S_2\right)\right)\right|\left(S_2-S_1\right)-F_{s1}$ (6)

$\begin{array}{c}{v}_2\frac{\text{d}{S}_2\left(t\right)}{\text{d}t}=\left|k_1\left(-{\alpha }_t\left(T_1-T_2\right)+{\beta }_s\left(S_1-S_2\right)\right)\right|\left(S_1-S_2\right)\\ +\left|k_2\left(-{\alpha }_t\left(T_3-T_2\right)+{\beta }_s\left(S_3-S_2\right)\right)\right|\left(S_3-S_2\right)+F_{s2}\end{array}$ (7)

$v_3\frac{\text{d}{S}_3\left(t\right)}{\text{d}t}=\left|k_2\left(-{\alpha }_t\left(T_3-T_2\right)+{\beta }_s\left(S_3-S_2\right)\right)\right|\left(S_2-S_3\right)-F_{s3}$ (8)

From (6) - (8), assuming that

$F_{s2}=F_{s1}+F_{s3}$

In this section, we take the division to be $v_1=v_2$ . According to [38], the volume of the global ocean is approximately 4 times the volume of Atlantic. Thus:

$v_1=v_2=\frac{v_3}{6}$

Denote

$y_1={\beta }_s\text{Δ}{S}_1,y_2={\beta }_s\text{Δ}{S}_2,x_1={\alpha }_t\text{Δ}{T}_1,x_2={\alpha }_t\text{Δ}{T}_2,$

We have

$\frac{\text{d}{y}_1}{\text{d}t}=-\frac{2k}{v_1}\left|y_1-x_1\right|y_1-\frac{k}{v_1}\left|y_2-x_2\right|y_2-\frac{\left(F_{s1}+F_{s2}\right){\beta }_s}{v_1}$ (9)

$\frac{\text{d}{y}_2}{\text{d}t}=-\frac{7k}{6v_1}\left|y_2-x_2\right|y_2-\frac{k}{v_1}\left|y_1-x_1\right|y_1-\left(\frac{F_{s3}}{6v_1}+\frac{F_{s2}}{v_1}\right){\beta }_s$ (10)

3.2. Steady States and Stability Analysis

In a steady state, we have

$\frac{\text{d}{y}_1}{\text{d}t}=0,$

$\frac{\text{d}{y}_2}{\text{d}t}=0,$

Noticing that $F_{s2}=F_{s1}+F_{s3}$ ,

From (9) − 2 * (10), we have

$0=\frac{\text{d}{y}_1}{\text{d}t}-2\frac{\text{d}{y}_2}{\text{d}t}=\frac{4}{3}\left(\frac{k}{v_1}\left|y_2-x_2\right|y_2+\frac{{\beta }_s\cdot F_{s3}}{v_1}\right)$ (11)

Similarly,

$\frac{k}{v_1}\left|y_1-x_1\right|y_1+\frac{F_{s1}\cdot {\beta }_s}{v_1}=0$ (12)

Thus, we have separated the entire large system into two independent subsystems. Recall that in the two-box model, there are three solutions when $x^2>4\frac{{\beta }_s{F}_s}{k}$ and one solution when $x^2<4\frac{{\beta }_s{F}_s}{k}$ . Similar to this approach, we categorize the problem into the following four cases.

Case 1: $x_1^2>4\frac{{\beta }_s{F}_{s1}}{k}$ and $x_2^2>4\frac{{\beta }_s{F}_{s2}}{k}$

Similar to the two-box case, denote:

$y_i^{\left(1\right)}=\frac{x}{2}-\frac{1}{2}\sqrt{x^2+4\frac{{\beta }_s{F}_{si}}{k}},$

$y_i^{\left(2\right)}=\frac{x}{2}-\frac{1}{2}\sqrt{x^2-4\frac{{\beta }_s{F}_{si}}{k}},$

$y_i^{\left(3\right)}=\frac{x}{2}+\frac{1}{2}\sqrt{x^2-4\frac{{\beta }_s{F}_{si}}{k}},$

for $i=1,2$ .

(We will also use this notation in Case 2, Case 3, and Case 4.)

It can be verified that $y_i^{\left(j\right)},j=1,2,3$ are the solutions for the function $\frac{k}{v_i}\left|y_i-x_i\right|y_i+\frac{F_{si}\cdot {\beta }_s}{v_i}=0,i=1,2$ .

As a result, there are up to $3\times 3=9$ steady states, i.e., combinations of $\left(y_1^{\left(j_1\right)},y_2^{\left(j_2\right)}\right)$ for $j_1,j_2=1,2,3$ .

The Jacobian analysis shows that 4 of these steady states are stable:

$\left(y_1,y_2\right)=\left(y_1^{\left(1\right)},y_2^{\left(1\right)}\right),\left(y_1^{\left(1\right)},y_2^{\left(3\right)}\right),\left(y_1^{\left(3\right)},y_2^{\left(1\right)}\right),\left(y_1^{\left(3\right)},y_2^{\left(3\right)}\right).$

And they represent the four main modes of the circulation, see Figure 5.

Case 2: $x_1^2<4\frac{{\beta }_s{F}_{s1}}{k}$ and $x_2^2>4\frac{{\beta }_s{F}_{s2}}{k}$

Similarly, it can be inferred that there are 3 steady states in total:

$\left(y_1,y_2\right)=\left(y_1^{\left(1\right)},y_2^{\left(j\right)}\right),j=1,2,3,$

with 2 of them being stable:

$\left(y_1,y_2\right)=\left(y_1^{\left(1\right)},y_2^{\left(1\right)}\right),\left(y_1^{\left(1\right)},y_2^{\left(3\right)}\right).$

Figure 5. The four main modes of circulation in Stommel’s three-box model.

Case 3: $x_1^2>4\frac{{\beta }_s{F}_{s1}}{k}$ and $x_2^2<4\frac{{\beta }_s{F}_{s2}}{k}$

Similarly, it can be inferred that there are 3 steady states in total:

$\left(y_1,y_2\right)=\left(y_1^{\left(j\right)},y_2^{\left(1\right)}\right),j=1,2,3,$

with 2 of them being stable:

$\left(y_1,y_2\right)=\left(y_1^{\left(1\right)},y_2^{\left(1\right)}\right),\left(y_1^{\left(3\right)},y_2^{\left(1\right)}\right).$

Case 4: $x_1^2<4\frac{{\beta }_s{F}_{s1}}{k}$ and $x_2^2<4\frac{{\beta }_s{F}_{s2}}{k}$

Similarly, it can be inferred that there exists only one steady state:

$\left(y_1,y_2\right)=\left(y_1^{\left(1\right)},y_2^{\left(1\right)}\right)$

And the steady state is stable.

Up to this point, we have discussed the steady-state solutions of this ordinary differential equation. However, for a given initial value, will the solution trajectory ultimately converge to one of these steady-state points? The following proof guarantees this.

Theorem 1 (Stability Analysis) There exists $a_1,a_2\in R$ such that

$\underset{t\to \infty }{\text{lim}}\left|y_1-a_1\right|=0$ and $\underset{t\to \infty }{\text{lim}}\left|y_2-a_2\right|=0$ .

Proof 3.1 We only prove the case $x_1^2\ge 4\frac{F_{s1}\cdot {\beta }_s}{V},x_2^2\ge 4\frac{F_{s2}\cdot {\beta }_s}{V}$ :

Denote

$\begin{array}{l}f\left(y\right)=-\frac{k}{v}\left|\left(y-x_1\right)\right|y-\frac{F_{s1}\cdot {\beta }_s}{v},\\ g\left(y\right)=-\frac{k}{v}\left|\left(y-x_2\right)\right|y-\frac{F_{s2}\cdot {\beta }_s}{v},\end{array}$

Thus

$\begin{array}{l}\frac{\text{d}{y}_1}{\text{d}t}=2f\left(y_1\right)+g\left(y_2\right),\\ \frac{\text{d}{y}_2}{\text{d}t}=f\left(y_1\right)+\frac{7}{6}g\left(y_2\right).\end{array}$

First consider the system $\frac{\text{d}y}{\text{d}t}=2f\left(y\right)$ .

Similar to the two-box model, because $x_1^2\ge 4\frac{F_{s1}\cdot {\beta }_s}{v}$ , the system has three steady states. Let $b_1<b_2$ be the stable point of the system:

$\frac{\text{d}y}{\text{d}t}=2f\left(y\right).$

$b_3$ is the unstable steady point. (Review this in the two-box case.)

Let $c_1<c_2$ be the stable point of the system:

$\frac{\text{d}y}{\text{d}t}=2g\left(y\right).$

$c_3$ is the unstable steady point. (Review this in the two-box case.)

Using these notations, noticing that in the three box model, when $y_1\le b_1$ , $f\left(y_1\right)\ge 0$ ; $b_1\le y_1\le b_3$ , $f\left(y_1\right)\le 0$ , so when $y_1\le b_3$ :

$\frac{\text{d}}{\text{d}t}\left|y_1-b_1\right|=-2\left|f\left(y_1\right)\right|+g\left(y_2\right)\mathrm{sgn}\left(y_1-b_1\right)\le -2\left|f\left(y_1\right)\right|+\left|g\left(y_2\right)\right|$ (13)

Similarly:

-If $y_1\ge b_3$ , then $\frac{\text{d}}{\text{d}t}\left|y_1-b_2\right|\le -2\left|f\left(y_1\right)\right|+\left|g\left(y_2\right)\right|$ (14)

-If $y_2\le c_3$ , then $\frac{\text{d}}{\text{d}t}\left|y_2-c_1\right|\le -2\left|f\left(y_1\right)\right|+\left|g\left(y_2\right)\right|$ (15)

-If $y_2\ge c_3$ , then $\frac{\text{d}}{\text{d}t}\left|y_2-c_2\right|\le -2\left|f\left(y_1\right)\right|+\left|g\left(y_2\right)\right|$ (16)

Let

$h\left(t\right)=\text{min}\left\{\left|y_1-b_1\right|+b_3-2b_2+b_1,\left|y_1-b_2\right|\right\}+\text{min}\left\{\left|y_2-c_1\right|+c_3-2c_2+c_1,\left|y_2-c_2\right|\right\}$

(We take the form $\text{min}\left\{\left|y_1-b_1\right|+b_3-2b_2+b_1,\left|y_1-b_2\right|\right\}$ to ensure the function strictly decreases.)

So, $h$ is continuous, and

from (13) and (14):

$\frac{\text{d}}{\text{d}t}\text{min}\left\{\left|y_1-b_1\right|+b_3-2b_2+b_1,\left|y_1-b_2\right|\right\}\le -2\left|f\left(y_1\right)\right|+\left|g\left(y_2\right)\right|$

from (15) and (16):

$\frac{\text{d}}{\text{d}t}\text{min}\left\{\left|y_2-c_1\right|+c_3-2c_2+c_1,\left|y_2-c_2\right|\right\}\le -\frac{7}{6}\left|g\left(y_2\right)\right|+\left|f\left(y_1\right)\right|$

Thus:

$\begin{array}{c}\frac{\text{d}h\left(t\right)}{\text{d}t}\le -\frac{1}{6}\left|g\left(y_2\right)\right|-\left|f\left(y_1\right)\right|\\ \le -\frac{1}{13}\left(\left|2f\left(y_1\right)+g\left(y_2\right)\right|+\left|f\left(y_1\right)+\frac{7}{6}g\left(y_2\right)\right|\right)\\ =-\frac{1}{13}\left(\left|\frac{\text{d}{y}_1}{\text{d}t}\right|+\left|\frac{\text{d}{y}_2}{\text{d}t}\right|\right)\end{array}$

Since $\frac{\text{d}h\left(t\right)}{\text{d}t}\le 0$ , and noticing that:

$h\left(t\right)\ge -\left|b_3-2b_2+b_1|-|c_3-2c_2+c_1\right|,$

we have:

${\displaystyle {\int }_0^{\infty }\left|\frac{\text{d}{y}_1}{\text{d}t}\right|\text{d}t}+{\displaystyle {\int }_0^{\infty }\left|\frac{\text{d}{y}_2}{\text{d}t}\right|\text{d}t}\le 13\left(\left|h\left(0\right)\right|+\left|b_3-2b_2+b_1\right|+\left|c_3-2c_2+c_1\right|\right)$

Therefore, there exists $a_1,a_2$ such that:

$\underset{t\to \infty }{\text{lim}}\left|y_1-a_1\right|=0$ and $\underset{t\to \infty }{\text{lim}}\left|y_2-a_2\right|=0.$

3.3. Discussion of the Extending Model

From the previous section, we can see that for a steady state, the number of distinct value q1 may take is completely determined by $F_{s1}$ .

We aim to study the solutions of $\text{Δ}{S}_1$ and $\text{Δ}{S}_2$ when given $F_{s1}$ , $F_{s3}$ , and initial values of $\text{Δ}{S}_1$ , $\text{Δ}{S}_2$ , as well as which steady-state solutions they will ultimately converge to. Here, we find that the influence of $F_{s3}$ on the entire system is non-negligible. Specifically, for given initial values and $F_{s1}$ , even when only $F_{s3}$ is varied, the final state of $\text{Δ}{S}_1$ may exhibit fundamentally different characteristics.

For instance, when setting $F_{s1}=1\text{m}/\text{y}$ with $F_{s3}=0.18\text{m}/\text{y}$ , compared to $F_{s1}=1\text{m}/\text{y}$ and $F_{s3}=0.2\text{m}/\text{y}$ , we notice great dissimilarity of the $F_{s1}$ temporal evolution curves, see Figure 6.

We notice that though when we take $F_{s1}=1,F_{s3}=0.2$ , AMOC collapse, contrary to that of $F_{s1}=1,F_{s3}=0.18$ .

Next, we show the minimum $F_{s1}$ value required to induce AMOC collapse under the initial conditions of $\text{Δ}{S}_3=-1$ , $\text{Δ}{S}_1$ at its first steady state, and $F_{s3}$ increasing from 0 to 0.3, see Figure 7.

We note that when $F_{s3}<0.186$ , the critical value of $F_{s1}$ that causes circulation 1 to collapse is 1.186. This occurs precisely because when $F_{s1}>1.186$ , in a steady state, circulation 1 has only one steady-state solution.

When we slowly increase the number of $F_{s3}$ , passing by 1.86, the corresponding $F_{s1}$ decreases sharply. Similar to the idea in [26], when we gradually increase $F_{s3}$ , the overall number of solutions switches from 3 to 1, so $\text{Δ}{S}_2$ should switch to that solution.

Recalling the equation $\frac{\text{d}{y}_1}{\text{d}t}=-\frac{2k}{v_1}\left|y_1-x_1\right|y_1-\frac{k}{v_2}\left|y_2-x_2\right|y_2-\frac{\left(F_{s1}+F_{s2}\right){\beta }_s}{v_1}$ . When $\text{Δ}{S}_2$ gradually decreases from its initial value to that solution, $\left(\frac{k}{v_1}\left|y_2-x_2\right|y_2+\frac{{\beta }_s\cdot F_{s3}}{v_1}\right)<0$ . Heuristically speaking, this negative term becomes the dominant part on right hand side, making $\frac{\text{d}{y}_1}{\text{d}t}<0$ , thereby driving $\text{Δ}{S}_1$ to decrease towards its second steady state point (see Figure 8 for the steady states of $q_1$ and $q_2$ corresponding to $F_{s1}$ and $F_{s3}$ ).

When $F_{s3}>1.86$ , the relationship is subtle. From observation, the graph of the function is almost a straight line. The least squares method was applied to approximate this plot (see Figure 9), from which we derived the expression

$F_{s1}=1.3655-1.9846F_{s3}.$

(a) $F_{s1}=1.00\text{m}/\text{y},F_{s3}=0.20\text{m}/\text{y}$ (b) $F_{s1}=1.00\text{m}/\text{y},F_{s3}=0.18\text{m}/\text{y}$

Figure 6. The projected trajectory of AMOC for $F_{s3}$ corresponding to 0.20 m/y and 0.18 m/y.

Figure 7. The the minimum value of $F_{s1}$ to make AMOC collapse as a function of $F_{s3}$ . (a) Given $F_{s3}$ , the minimum value of $F_{s1}$ to make AMOC collapse; (b) A 3-D perspective showing the final strength of q1(AMOC) for different$F_{s1},F_{s3}$ .

(a) The steady states of $q_1$ corresponding to $F_{s1}$ (b) The steady states of $q_2$ corresponding to $F_{s3}$

Figure 8. Steady states of salinity difference and AMOC as a function of $F_s$ .

Figure 9. The least square fit.

As we have seen that the introduction of a third box enables the AMOC to “bypast” one of the robust steady states present in the two-box model. The physical extrapolation of this result reveals that the Pacific becomes a dynamic booster, imparting sufficient inertia to the system so that it cannot settle into the original potential wells.

In the traditional two-box model, transitioning from one steady state to another requires overcoming an energy barrier, acting as a kind of “sticky” hysteresis. However, in the three-box model, the inclusion of the Pacific introduces a huge volume and salinity capacity, endowing the system with a “kinetic flywheel” effect. When the system begins to evolve from one state to another under external forcing, the established salinity gradient between the Atlantic and Pacific drives sustained circulation and exchange. The momentum of this process prevents the system from pausing at the intermediate equilibrium points that existed before—much like a swing being pushed continuously does not come to rest at the midpoint but swings directly from one high point to another, thereby “bypass” the seemingly robust intermediate steady states that would have existed without the Pacific.

This phenomenon can be analogized to the formation process of “explosive cyclogenesis” in weather systems. During cyclogenesis, when certain conditions reach a critical threshold, the system does not slowly approach equilibrium. Instead, it gains immense inertia, rapidly skipping over weak stable phases and entering a period of rapid intensification or “explosion”. Similarly, in our three-box model, the additional capacity provided by the third box provides the AMOC with similar inertia during its state evolution, causing it to bypass its original steady-state structure.

4. A Simple Model for Dansgaard-Oeschger Cycles

4.1. Construction of the Model

As is well known from modelling, the Atlantic Meridional Overturning Circulation (AMOC) has more than one steady-state mechanism available. During D-O events, the vigorous AMOC brings increased quantities of warm water mass into the North Atlantic Ocean and Arctic Ocean, increasing the freshwater transport from Arctic to Northern Hemisphere landmasses as well as the subsequent, amplified freshwater forcing [39] [40]. As discussed in the Stommel’s Model in Section 2, when the amplitude of the freshwater forcing passes some critical threshold, the stable, strong, steady state vanishes, and the AMOC switches to the weak steady state. In this weak state, the northward heat transport is substantially reduced relative to the strong state, leading to a decrease in freshwater runoff and hence freshwater forcing. If the forcing becomes weak enough for the upstream process to cease and fresh water no longer moves northward, the strong AMOC state can be reinitiated, bringing about a complete cyclical transition through this means between the two distinct states.

In recent research [36], the authors reconstructed a two-box AMOC model by introducing a temporal delay in the oceanic transport between the boxes through a delay function, demonstrating the self-sustained multidecadal variability of AMOC. We hypothesize that the mechanism of Dansgaard-Oeschger (D-O) events may share similar principles. Notably, the polar see-saw hypothesis discussed in [12] [14] suggests that during deepwater formation in polar regions, substantial changes (e.g., during intense glacial growth or depletion) require prolonged periods to manifest their effects in the opposite hemisphere.

This implies that when a strong AMOC resumes, Arctic runoff reduction shows delayed response. Consequently, the strong AMOC persists while continually “charging” the Arctic ice system with energy. When the AMOC’s influence eventually reaches the Arctic, the resultant dramatic increase in runoff forces the AMOC to collapse into a weak state. Due to this accumulated energy reserve, the elevated Arctic runoff persists even after AMOC weakening, maintaining the weak AMOC phase until the stored energy is fully depleted.

We suppose the volume of the volume, salinity and temperature for the Arctic are $V_1,S_1,T_1$ , respectively, and those of the high latitude Atlantic and low latitude Atlantic are $V_2,S_2,T_2$ and $V_3,S_3,T_3$ . Given the extended duration of Dansgaard-Oeschger (D-O) events, we explicitly consider both temperature and salinity as dynamic variables in our model framework. See Figure 10 for our model.

Figure 10. Schematic diagram of the model.

Our analysis begins with the freshwater forcing. The freshwater forcing primarily consists of two components: runoff from ice sheets and net rainfall. We assume the fresh water force for the Arctic and low-latitude Atlantic is dominated by the runoff of the ice-sheet and net rainfall, respectively, while the freshwater force for the high-latitude is dominated by both runoff from ice-sheet and net rainfall.

Firstly, let’s assess the effect of ice-sheet runoff on marine systems. Let R represent the runoff rate, and given freshwater’s near-zero salinity and temperature, we establish the core salinity and heat balance equation for a simple box:

$v\frac{\text{d}S}{\text{d}t}=-SR,v\frac{\text{d}T}{\text{d}t}=-TR,$

In reality, since the runoff R is much smaller than the ocean volumes, we assume that ocean volumes are constants. The volume of the Arctic Ocean is approximately 1.9 × 10⁷ km³. The volume of the North Atlantic Ocean is approximately 1.6 × 10⁸ km³ [38]. In this section, we define the partition as $v_1=4\times {10}^7{\text{km}}^3$ and $v_2=1.2\times {10}^8{\text{km}}^3$ .

Besides, in D-O event, fresh water runs to the high latitude North Atlantic and the Arctic [41], denote them by $R_N,R_A$ , respectively. We take the ratio to be 2:1:

$R_A=\frac{R}{3},R_N=\frac{2}{3}R$

Again, let’s consider the effect of net rainfall. For simplicity, we assume the effect of net rainfall on high-latitude North Atlantic and low-latitude North Atlantic is constants. Denote it by $-Fs$ and $Fs$ , respectively.

Let’s move on to the effect of AMOC on the ice-sheet runoff R. Varying meltwater runoff to the North Atlantic can occur in association with fluctuations of the ice margin [33]. There are a number of inter-related considerations here. We will use the result in [33] that the runoff R will depend linearly on AMOC. Besides, our framework incorporates a temporal lag (denoted as $d_1$ ) between the low-latitude Atlantic meridional overturning circulation and its impact on ice-sheet runoff dynamics:

$R=R_0+a\cdot q_2\left(t-d_1\right)$

Estimates of the period of temporal lag vary and one typical estimate is 400 years and we will illustrate it in the next section. Regarding the selection of parameters $R_0,a$ and $Q_0$ , we refer to the work in [35]. We will use $R_0~0.1\text{Sv}$ , $a~0.01$ and $Q_0~20\text{Sv}$ .

For the volume transport $q_1,q_2$ , we assume that they are proportional to the density difference between the boxes. Additionally, we incorporate another temporal lag (denoted as $d_2$ ) to characterize the mean advective time delay, a choice for $d_2$ is about 20 years:

$q_1=k_1\left(-{\alpha }_t\left(T_1\left(t-d_2\right)-T_2\left(t-d_2\right)\right)+{\beta }_s\left(S_1\left(t-d_2\right)-S_2\left(t-d_2\right)\right)\right)$

$q_2=k_2\left(-{\alpha }_t\left(T_2\left(t-d_2\right)-T_3\left(t-d_2\right)\right)+{\beta }_s\left(S_2\left(t-d_2\right)-S_3\left(t-d_2\right)\right)\right)$

The dynamic equations of salinity for our three-box model become:

$v_1\frac{\text{d}{S}_1\left(t\right)}{\text{d}t}=\left|q_1\right|\cdot \left(S_2\left(t-d_2\right)-S_1\left(t-d_2\right)\right)-S_1\left(t\right)\cdot R_A$

$\begin{array}{c}{v}_2\frac{\text{d}{S}_2\left(t\right)}{\text{d}t}=\left|q_1\right|\cdot \left(S_1\left(t-d_2\right)-S_2\left(t-d_2\right)\right)\\ +\left|q_2\right|\cdot \left(S_3\left(t-d_2\right)-S_2\left(t-d_2\right)\right)-S_2\left(t\right)\cdot R_N-F_s\end{array}$

$v_3\frac{\text{d}{S}_3\left(t\right)}{\text{d}t}=\left|q_2\right|\cdot \left(S_2\left(t-d_2\right)-S_3\left(t-d_2\right)\right)+F_s$

In addition, we consider temperature as dynamic variables in our model framework. Assuming the ambient temperature for Arctic, high-latitude North Atlantic and low-latitude North Atlantic are $c_1,c_2,c_3$ , and let $r_1,r_2,r_3$ characterizing the speed of heat conducting into the water, we derive:

$v_1\frac{\text{d}{T}_1\left(t\right)}{\text{d}t}=\left|q_1\right|\cdot \left(T_2\left(t-d_2\right)-T_1\left(t-d_2\right)\right)-T_1\left(t\right)\cdot R_A+r_1\left(c_1-T_1\left(t\right)\right)$

$\begin{array}{c}{v}_2\frac{\text{d}{T}_2\left(t\right)}{\text{d}t}=\left|q_1\right|\cdot \left(T_1\left(t-d_2\right)-T_2\left(t-d_2\right)\right)+\left|q_2\right|\cdot \left(T_3\left(t-d_2\right)-T_2\left(t-d_2\right)\right)\\ -T_2\left(t\right)\cdot R_N+r_2\left(c_2-T_2\left(t\right)\right)\end{array}$

$v_3\frac{\text{d}{T}_3\left(t\right)}{\text{d}t}=\left|q_2\right|\cdot \left(T_2\left(t-d_2\right)-T_3\left(t-d_2\right)\right)+r_3\left(c_3-T_3\left(t\right)\right)$

4.2. Explanation of Temporal Lag

The delayed feedback mechanisms within the Arctic climate system exhibit unique temporal characteristics during abrupt climate transitions. When the AMOC intensifies, warm water reaches the Arctic marine margins within 1 - 3 years [42]. This rapid advection stems from the Last Glacial Maximum paleogeographic configuration: lower sea levels extended ice sheets to continental shelf edges, enabling direct thermal erosion of ice-shelf basal zones—a stark contrast to modern fjord-dominated retreat [43].

Arctic ice sheets respond to oceanic forcing through three distinct stages:

1) Initial basal melting (0 - 15 years): Subsurface warming triggers ice-shelf thinning at rates > 10 m/y, recorded by abrupt increases in ice-rafted debris with Laurentide-sourced dolomite in North Atlantic sediments [44];

2) Ice dynamic acceleration (15 - 30 years): Reduced buttressing causes inland ice-stream acceleration, evidenced by layer-thinning in Greenland ice cores [45];

3) Runoff peak (30 - 50 years): Combined meltwater discharge and iceberg calving release rates exceeded 0.1 Sv freshwater, as reported in [46] [47].

In addition, the accumulation of freshwater is subject to complex spatiotemporal patterns, with proglacial lakes attaining outburst capacities (up to 500 Gt) in 50 - 150 years, as seen in varved sediment records from Lake Agassiz [48], and subglacial aquifers storing freshwater on time scales of 200 to 500 years, as evidenced by porewater Cl⁻ profiles and ³⁶Cl dating [49].

Despite decadal-scale freshwater discharge into the ocean, substantial weakening of the AMOC manifests after approximately 250 - 350 years. This delayed response primarily stems from time-dependent accumulation processes in critical convection zones. Norwegian Sea sediment Nd records reveal a 260 ± 40 yr delay between ice-rafted debris events and subsequent AMOC collapse [50].

There is a high-resolution archive of documents in lag hierarchy as follows. Antarctic EPICA ice core δ18O excursions lag Greenland stadial—interstadial transitions by 310 ± 50 years on average [51]. South Atlantic Nd data further indicate that AMOC signal propagation via oceanic pathways after approximately 190 ± 30 years, giving rise to a bipolar see-saw lag of 300 ± 70 years [52].

By integrating all feedback stages into one sequence of assessment, we suggest a delay of 400 years as a very conservative estimate representing a robust time-lag to ice-sheet-ocean-climate interaction throughout the system.

4.3. The Numerical Results of Our Model

Self-sustained multidecadal AMOC oscillations are obtained in the revised three-box Stommel model. Figure 11 shows the numerical result at $d_1=400y,d_2=20y$ .

Figure 11. Numerical solutions of our model. (a) The temperature of the Arctic, high-latitude Atlantic and low-latitude Atlantic over the years 0 - 20 kyr; (b) The strength of AMOC and circulation between the high-latitude box and the Arctic over the time series years 0 - 20 kyr.

4.4. Self-Sustained Oscillations

First, let’s examine the temperature variations. When the strong AMOC initiates, the high-latitude North Atlantic and Arctic regions experience abrupt warming, with the former showing a rapid temperature increase of approximately 6˚C and the latter rising by about 4˚C within a short timeframe, which is consistent with the results shown in [35]. Following the collapse of the strong AMOC, temperatures in both regions gradually decrease and eventually return to their initial levels.

And from the graphs, we can observe that a Dansgaard-Oeschger (D-O) cycle consists of the following phases:

1) At point a (see Figure 12), as freshwater forcing ($R_N,R_A$ ) declines to a critical threshold, it triggers the activation of a strong AMOC.

2) Point b marks the peak intensity of the strong AMOC. After about 400 years, under persistent ice sheet runoff from Arctic ice sheets (primarily meltwater input), freshwater forcing increases, leading to the collapse of the strong AMOC and its transition to a weak mode. Point c denotes the terminus of the strong AMOC (Strong phase).

3) From point d to e, under weak AMOC conditions, ice sheet runoff diminishes, and freshwater forcing gradually decreases. During this period, freshwater forcing declines to negative values, leaving the system in a monostable state dominated solely by the weak AMOC. Notably, the AMOC intensity increases with periodic fluctuations. This pattern primarily arises because when Arctic runoff (R) remains stable, the AMOC can reach its equilibrium state within a timescale that is a bit shorter than the characteristic temporal lag of the system. The cumulative effects of AMOC’s influence on Arctic runoff only become manifest after the subsequent 400-year phase. This delayed feedback triggers freshwater forcing adjustments, driving the AMOC toward a new equilibrium state.

4) At point e, freshwater forcing declines to negative value which triggers the next cycle.

Figure 12. Modelled cycle about the AMOC.

All the cycles go through the same sequence and each lasts around 1500 - 2000 years, which is about 4 - 5 times the typical time delay of the entire system. The cycle duration largely relies on Phase C—the most time-consuming phase—lasting for a very long period of time.

In our model, the influence of AMOC (Atlantic Meridional Overturning Circulation) on ice-sheet runoff only becomes evident after 400 years. This temporal lag allows the AMOC to reach its new steady state under the ice-sheet runoff values from the previous phase, consequently causing the AMOC variation to demonstrate this pattern. Therefore, we can approximately discretize the model

Figure 13. The discretized model (one step).

regarding time t into a recursive sequence(see Figure 13). Strong AMOC can only restart when freshwater forcing drops below zero. Consequently, the timing of strong AMOC restart depends on which term in this sequence first attains negative freshwater forcing values. Figure 13 illustrates this process.

5. Conclusions

We have developed a three-box Stommel model to simulate the mechanisms of the Atlantic Meridional Overturning Circulation (AMOC). The principle is straightforward: the global ocean is partitioned into three adjacent boxes—the high-latitude North Atlantic, the low-latitude North Atlantic, and the rest of the global ocean. Temperature and salinity gradients between these regions drive the circulation. Our model relies on key assumptions: 1) the freshwater forcing applied to each box is constant, and 2) the temperature differences between regions remain fixed. Future work may involve relaxing these assumptions for deeper exploration.

Our findings reveal that the third box (representative of the global ocean) plays an important part in steering AMOC to specific steady-state solutions; if it receives stronger freshwater forcing, it tends to bypass robust steady states, and gradually converge to weak equilibria which never happens with the conventional two-box model and suggests that the AMOC collapse under the RCP 8.5 warming scenario would take place much earlier than thought before.

Subsequently, based on the coupled time-delayed feedback of the three-box Stommel equations, we constructed a D-O cycles model. The two main factors are: 1) Bistability Loop—Strong AMOC enhances the ice sheet’s freshwater release and eventually leads to the collapse of the ice sheet, whereas weak AMOC weakens the freshwater forcing and promotes the restart of the AMOC. 2) Delayed response: A temporal lag (400 yr) exists between AMOC formation and its impact on ice-sheet mass balance.

This framework effectively captures key characteristics of observed Dansgaard-Oeschger events, including: Self-sustained oscillations between AMOC states, asymmetric timescale between cooling/warming phases and abrupt warming amplitudes. This motivates us to consider whether similar delayed response mechanisms exist when investigating analogous self-sustained oscillation events in future studies.

Our model development aims to overcome three limitations: one is that it must incorporate stochastic forcing to represent the effect of climate noise [53], the second is that a better representation of NADW formation zones requires higher spatial resolution [54], and the third is to adopt full-coupled ice sheet dynamics as in Asay-Davis et al. [55]. Further improved physics with retained delayed feedback paradigm should be expected after we implement those changes.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Broeker, W.S. (1991) The Great Ocean Conveyor. Oceanography, 4, 79-89.[CrossRef]
[2]	Broecker, W.S., Peteet, D.M. and Rind, D. (1985) Does the Ocean-Atmosphere System Have More than One Stable Mode of Operation? Nature, 315, 21-26.[CrossRef]
[3]	Manabe, S. and Stouffer, R.J. (1997) Coupled Ocean-Atmosphere Model Response to Freshwater Input: Comparison to Younger Dryas Event. Paleoceanography, 12, 321-336.[CrossRef]
[4]	Rahmstorf, S. (2002) Ocean Circulation and Climate during the Past 120,000 Years. Nature, 419, 207-214.[CrossRef] [PubMed]
[5]	Zhang, R., Sutton, R., Danabasoglu, G., Kwon, Y., Marsh, R., Yeager, S.G., et al. (2019) A Review of the Role of the Atlantic Meridional Overturning Circulation in Atlantic Multidecadal Variability and Associated Climate Impacts. Reviews of Geophysics, 57, 316-375.[CrossRef]
[6]	Galochkina, M., Makarova, M., Miller, K.G., Browning, J.V., Keating, R.S. and Wright, J.D. (2023) Multispecies Planktonic and Benthic Foraminiferal Stable Isotopes from North Atlantic Subtropical Site 558: Thermocline Intensification during the Mid-Miocene Climate Transition. Journal of Foraminiferal Research, 53, 143-156.[CrossRef]
[7]	Andersen, K.K. (2004) High-Resolution Record of Northern Hemisphere Climate Extending into the Last Interglacial Period. Nature, 431, 147-151.[CrossRef] [PubMed]
[8]	Dansgaard, W., Johnsen, S.J., Clausen, H.B., Dahl-Jensen, D., Gundestrup, N.S., Hammer, C.U., et al. (1993) Evidence for General Instability of Past Climate from a 250-Kyr Ice-Core Record. Nature, 364, 218-220.[CrossRef]
[9]	Johnsen, S.J., Clausen, H.B., Dansgaard, W., Fuhrer, K., Gundestrup, N., Hammer, C.U., et al. (1992) Irregular Glacial Interstadials Recorded in a New Greenland Ice Core. Nature, 359, 311-313.[CrossRef]
[10]	Kelsey, A.M., et al. (2015) An Astronomical Correspondence to the 1470 Year Cycle of Abrupt Climate Change. Climate of the Past Discussions.
[11]	Wolff, E.W., Chappellaz, J., Blunier, T., Rasmussen, S.O. and Svensson, A. (2010) Millennial-Scale Variability during the Last Glacial: The Ice Core Record. Quaternary Science Reviews, 29, 2828-2838.[CrossRef]
[12]	Barker, S., Knorr, G., Edwards, R.L., Parrenin, F., Putnam, A.E., Skinner, L.C., et al. (2011) 800,000 Years of Abrupt Climate Variability. Science, 334, 347-351.[CrossRef] [PubMed]
[13]	Stocker, T.F. and Johnsen, S.J. (2003) A Minimum Thermodynamic Model for the Bipolar Seesaw. Paleoceanography, 18, Article No. 1087.[CrossRef]
[14]	Severinghaus, J.P. (2009) Southern See-Saw Seen. Nature, 457, 1093-1094.[CrossRef] [PubMed]
[15]	Bryden, H.L., Longworth, H.R. and Cunningham, S.A. (2005) Slowing of the Atlantic Meridional Overturning Circulation at 25˚ N. Nature, 438, 655-657.[CrossRef] [PubMed]
[16]	Caesar, L., McCarthy, G.D., Thornalley, D.J.R., Cahill, N. and Rahmstorf, S. (2021) Current Atlantic Meridional Overturning Circulation Weakest in Last Millennium. Nature Geoscience, 14, 118-120.[CrossRef]
[17]	Rahmstorf, S., Box, J.E., Feulner, G., Mann, M.E., Robinson, A., Rutherford, S., et al. (2015) Exceptional Twentieth-Century Slowdown in Atlantic Ocean Overturning Circulation. Nature Climate Change, 5, 475-480.[CrossRef]
[18]	Smeed, D.A., Josey, S.A., Beaulieu, C., Johns, W.E., Moat, B.I., Frajka-Williams, E., et al. (2018) The North Atlantic Ocean Is in a State of Reduced Overturning. Geophysical Research Letters, 45, 1527-1533.[CrossRef]
[19]	Boers, N. (2021) Publisher Correction: Observation-Based Early-Warning Signals for a Collapse of the Atlantic Meridional Overturning Circulation. Nature Climate Change, 11, Article No. 1001.[CrossRef]
[20]	Ditlevsen, P. and Ditlevsen, S. (2023) Warning of a Forthcoming Collapse of the Atlantic Meridional Overturning Circulation. Nature Communications, 14, Article No. 4254.[CrossRef] [PubMed]
[21]	Ganopolski, A. and Rahmstorf, S. (2001) Rapid Changes of Glacial Climate Simulated in a Coupled Climate Model. Nature, 409, 153-158.[CrossRef] [PubMed]
[22]	Hawkins, E., Smith, R.S., Allison, L.C., Gregory, J.M., Woollings, T.J., Pohlmann, H., et al. (2011) Correction to “Bistability of the Atlantic Overturning Circulation in a Global Climate Model and Links to Ocean Freshwater Transport”. Geophysical Research Letters, 38, L16699.[CrossRef]
[23]	Johnson, H.L., Marshall, D.P. and Sproson, D.A.J. (2007) Reconciling Theories of a Mechanically Driven Meridional Overturning Circulation with Thermohaline Forcing and Multiple Equilibria. Climate Dynamics, 29, 821-836.[CrossRef]
[24]	Rahmstorf, S., Crucifix, M., Ganopolski, A., Goosse, H., Kamenkovich, I., Knutti, R., et al. (2005) Thermohaline Circulation Hysteresis: A Model Intercomparison. Geophysical Research Letters, 32, L23605.[CrossRef]
[25]	Shin, Y., et al. (2025) Reconciled Warning Signals in Observations and Models Imply Approaching AMOC Tipping Point. arXiv: 2503.22111.
[26]	Tziperman, E. (2022) Global Warming Science. Princeton University Press, 99-120.
[27]	Zhang, R., Follows, M. and Marshall, J. (2002) Mechanisms of Thermohaline Mode Switching with Application to Warm Equable Climates. Journal of Climate, 15, 2056-2072.[CrossRef]
[28]	Colin de Verdière, A. (2007) A Simple Model of Millennial Oscillations of the Thermohaline Circulation. Journal of Physical Oceanography, 37, 1142-1155.[CrossRef]
[29]	Colin de Verdière, A., Ben Jelloul, M. and Sévellec, F. (2006) Bifurcation Structure of Thermohaline Millennial Oscillations. Journal of Climate, 19, 5777-5795.[CrossRef]
[30]	Menviel, L., Timmermann, A., Friedrich, T. and England, M.H. (2014) Hindcasting the Continuum of Dansgaard-Oeschger Variability: Mechanisms, Patterns and Timing. Climate of the Past, 10, 63-77.[CrossRef]
[31]	Birchfield, G.E. and Broecker, W.S. (1990) A Salt Oscillator in the Glacial Atlantic? 2. A “Scale Analysis” Model. Paleoceanography, 5, 835-843.[CrossRef]
[32]	Broecker, W.S., Bond, G., Klas, M., Bonani, G. and Wolfli, W. (1990) A Salt Oscillator in the Glacial Atlantic? 1. The Concept. Paleoceanography, 5, 469-477.[CrossRef]
[33]	Clark, P.U., Marshall, S.J., Clarke, G.K.C., Hostetler, S.W., Licciardi, J.M. and Teller, J.T. (2001) Freshwater Forcing of Abrupt Climate Change during the Last Glaciation. Science, 293, 283-287.[CrossRef] [PubMed]
[34]	Zhang, X., Lohmann, G., Knorr, G. and Purcell, C. (2014) Abrupt Glacial Climate Shifts Controlled by Ice Sheet Changes. Nature, 512, 290-294.[CrossRef] [PubMed]
[35]	Wolff, E., et al. (2014) A Freshwater Starvation Mechanism for Dansgaard-Oeschger Cycles. 2014 AGU Fall Meeting, San Francisco,1 December 2014, PP11E-04.
[36]	Wei, X. and Zhang, R. (2022) A Simple Conceptual Model for the Self-Sustained Multidecadal AMOC Variability. Geophysical Research Letters, 49, e2022GL099800.[CrossRef]
[37]	Stommel, H. (1961) Thermohaline Convection with Two Stable Regimes of Flow. Tellus A: Dynamic Meteorology and Oceanography, 13, 224-230.[CrossRef]
[38]	Eakins, B.W. and Sharman, G.F. (2010) Volumes of the Worlds Oceans from ETOPO1. NOAA National Geophysical Data Center.
[39]	Jungclaus, J.H., Haak, H., Latif, M. and Mikolajewicz, U. (2005) Arctic-North Atlantic Interactions and Multidecadal Variability of the Meridional Overturning Circulation. Journal of Climate, 18, 4013-4031.[CrossRef]
[40]	Li, D., Zhang, R. and Knutson, T. (2018) Comparison of Mechanisms for Low-Frequency Variability of Summer Arctic Sea Ice in Three Coupled Models. Journal of Climate, 31, 1205-1226.[CrossRef]
[41]	Rignot, E., Mouginot, J., Scheuchl, B., van den Broeke, M., van Wessem, M.J. and Morlighem, M. (2019) Four Decades of Antarctic Ice Sheet Mass Balance from 1979-2017. Proceedings of the National Academy of Sciences, 116, 1095-1103.[CrossRef] [PubMed]
[42]	Zhang, R. and Vallis, G.K. (2007) The Role of Bottom Vortex Stretching on the Path of the North Atlantic Western Boundary Current and on the Northern Recirculation Gyre. Journal of Physical Oceanography, 37, 2053-2080.[CrossRef]
[43]	Bassis, J.N., Petersen, S.V. and Mac Cathles, L. (2017) Heinrich Events Triggered by Ocean Forcing and Modulated by Isostatic Adjustment. Nature, 542, 332-334.[CrossRef] [PubMed]
[44]	Rignot, E., Xu, Y., Menemenlis, D., Mouginot, J., Scheuchl, B., Li, X., et al. (2016) Modeling of Ocean-Induced Ice Melt Rates of Five West Greenland Glaciers over the Past Two Decades. Geophysical Research Letters, 43, 6374-6382.[CrossRef]
[45]	Alvarez-Solas, J., Charbit, S., Ritz, C., Paillard, D., Ramstein, G. and Dumas, C. (2010) Links between Ocean Temperature and Iceberg Discharge during Heinrich Events. Nature Geoscience, 3, 122-126.[CrossRef]
[46]	van den Broeke, M.R., Enderlin, E.M., Howat, I.M., Kuipers Munneke, P., Noël, B.P.Y., van de Berg, W.J., et al. (2016) On the Recent Contribution of the Greenland Ice Sheet to Sea Level Change. The Cryosphere, 10, 1933-1946.[CrossRef]
[47]	Goelzer, H., Nowicki, S., Payne, A., Larour, E., Seroussi, H., Lipscomb, W.H., et al. (2020) The Future Sea-Level Contribution of the Greenland Ice Sheet: A Multi-Model Ensemble Study of ISMIP6. The Cryosphere, 14, 3071-3096.[CrossRef]
[48]	Teller, J.T., Leverington, D.W. and Mann, J.D. (2002) Freshwater Outbursts to the Oceans from Glacial Lake Agassiz and Their Role in Climate Change during the Last Deglaciation. Quaternary Science Reviews, 21, 879-887.[CrossRef]
[49]	Person, M., Dugan, B., Swenson, J.B., Urbano, L., Stott, C., Taylor, J., et al. (2003) Pleistocene Hydrogeology of the Atlantic Continental Shelf, New England. Geological Society of America Bulletin, 115, 1324-1343.[CrossRef]
[50]	Galaasen, E.V., Ninnemann, U.S., Kessler, A., Irvalı, N., Rosenthal, Y., Tjiputra, J., et al. (2020) Interglacial Instability of North Atlantic Deep Water Ventilation. Science, 367, 1485-1489.[CrossRef] [PubMed]
[51]	Barker, S., Chen, J., Gong, X., Jonkers, L., Knorr, G. and Thornalley, D. (2015) Icebergs Not the Trigger for North Atlantic Cold Events. Nature, 520, 333-336.[CrossRef] [PubMed]
[52]	Skinner, L.C., Primeau, F., Freeman, E., de la Fuente, M., Goodwin, P.A., Gottschalk, J., et al. (2017) Radiocarbon Constraints on the Glacial Ocean Circulation and Its Impact on Atmospheric CO2. Nature Communications, 8, Article No. 16010.[CrossRef] [PubMed]
[53]	Ganopolski, A. and Rahmstorf, S. (2001) Rapid Changes of Glacial Climate Simulated in a Coupled Climate Model. Nature, 409, 153-158.[CrossRef] [PubMed]
[54]	Zhang, X., Lohmann, G., Knorr, G. and Xu, X. (2013) Different Ocean States and Transient Characteristics in Last Glacial Maximum Simulations and Implications for Deglaciation. Climate of the Past, 9, 2319-2333.[CrossRef]
[55]	Asay-Davis, X.S., Cornford, S.L., Durand, G., Galton-Fenzi, B.K., Gladstone, R.M., Gudmundsson, G.H., et al. (2016) Experimental Design for Three Interrelated Marine Ice Sheet and Ocean Model Intercomparison Projects: MISMIP V. 3 (MISMIP +), ISOMIP V. 2 (ISOMIP +) and MISOMIP V. 1 (misomip1). Geoscientific Model Development, 9, 2471-2497. [Google Scholar] [CrossRef]