2 Overview of AMOC and the Stommel 2-box model

2.1 The Stommel 2-box model

A well-known class of climate models, where salinity-dominated and temperature-dominated states are bi-stable, is that of thermohaline circulation (THC). Here, abrupt qualitative changes are possible, see Alley et al. [Reference Alley2], Marotzke [Reference Marotzke15] or Rahmstorf [Reference Rahmstorf18, Reference Rahmstorf19]. Recently, Rahmstorf was able to find evidence of weakening occurring around these abrupt changes in the AMOC system of ocean patterns [Reference Caesar, Rahmstorf, Robinson, Feulner and Saba7]. This evidence of ocean dynamics responding to changes in surface temperature underscores the need to understand the transitions in these types of systems. We note that such transitions can be either smooth or nonsmooth (as described in [Reference Bernardo, Budd, Champneys and Kowalczyk4]). In this paper, we focus on the commonly used Stommel 2-box model [Reference Stommel23] as an exemplar for studying the transitions in the THC (or more generally, the dynamical impact of NSF bifurcations between equilibrium states) in a realistic climate model. We begin with the nondimensionalized Stommel 2-box model as given in [Reference Dijkstra8]:

(2.1) $$ \begin{align} \begin{aligned} \dot{{\cal T}} & = \eta_1-{\cal T}(1+|{\cal T}-{\cal S}|), \\ \dot{{\cal S}} & = \eta_2-{\cal S}(\eta_3+|{\cal T}-{\cal S}|). \end{aligned} \end{align} $$

Here, the variables ${\cal T}$ and ${\cal S}$ are the dimensionless equatorial-to-pole differences for temperature and salinity, respectively. The parameters $\eta _1$ , $\eta _2$ and $\eta _3$ are also dimensionless quantities, with $\eta _1$ representing thermal variation, $\eta _2$ is the freshwater flux, and $\eta _3$ is the ratio of relaxation times of temperature and salinity. The dimensionless AMOC strength is captured by the difference

$$ \begin{align*}{\cal V}= {\cal T}-{\cal S},\end{align*} $$

which plays an important role throughout the dynamical analysis. We can then express this system as

(2.2) $$ \begin{align} \begin{aligned} \dot{{\cal T}} & = \eta_1-{\cal T}(1+|{\cal V}|), \\ \dot{{\cal V}} & = \eta_1 -\eta_2 +\eta_3({\cal T}-{\cal V})-{\cal T}-{\cal V}|{\cal V}|. \end{aligned} \end{align} $$

With the dependence on the absolute value $|{\cal V}|$ , (2.2) is a nonsmooth dynamical system. It has a discontinuity surface  $\Sigma $ given by

$$ \begin{align*} \Sigma = \{{\cal V} =0\}, \end{align*} $$

across which we see a discontinuity in $d^2{\cal T}/dt^2$ and $d^2{\cal V}/dt^2.$ The equations for ${\cal T}$ and ${\cal V}$ then describe different dynamics in $\Sigma ^+$ (the temperature-dominated state) and $\Sigma ^-$ (the salinity-dominated state) for

$$ \begin{align*} \Sigma^+ = \{({\cal T},{\cal V})\mid {\cal V}> 0 \}, \quad \Sigma^- = \{({\cal T},{\cal V})\mid {\cal V}< 0\}. \end{align*} $$

The model is nonsmooth through the action of the nonlinearity $|{\cal T} - {\cal S}|$ and takes the form of a piecewise-smooth continuous system with a degree of discontinuity equal to 2 [Reference Bernardo, Budd, Champneys and Kowalczyk4, Reference Simpson22].

2.2 Static dynamics

It is straightforward to analyse the Stommel 2-box model in the case of static parameters $\eta _i$ [Reference Dijkstra8, Reference Kaper and Engler12]. A standard analysis of the static model, where typically $\eta _1$ and $\eta _3$ are fixed, and $\eta _2$ is treated as a bifurcation parameter, yields stability regions for the temperature and salinity-dominated states. Taking values of $\eta _1$ and $\eta _3$ (typically $\eta _1 = 3, \eta _3 = 0.3$ ) as is usual in applications [Reference Dijkstra8], there are either three or one fixed points. In the case of three fixed points, we identify two different critical bifurcation points, denoted as $\eta _{2sn}$ and $\eta _{2sf} = \eta _1 \eta _3$ , with $\eta _{2sn}>\eta _{2sf}$ .

For $\eta _{2sn}>\eta _2>\eta _{2sf}$ , there are two fixed points in $\Sigma ^+$ which are respectively a saddle (S) and a stable node (N). The saddle and node coalesce at the (smooth) SNB $\eta _{2sn}$ and cease to exist for $\eta _2> \eta _{2sn}$ .

For $\eta _2>\eta _{2sf}$ , there is a further fixed point in $\Sigma ^-$ , which is a stable focus (F). If $\eta _2< \eta _{2sf}$ , there is only the stable node N in $\Sigma ^+$ . These observations are illustrated in Figure 1 for ${\cal V}$ versus $\eta _2$ . Note that N corresponds to the temperature-dominated state with ${\cal T}> {\cal S}$ , and F corresponds to the salinity-dominated state with ${\cal S}> {\cal T}$ .

Figure 1 The static bifurcation diagram for the Stommel 2-box model showing ${\cal V}$ as a function of $\eta _2$ (blue) when $\eta _1 = 3, \eta _3 = 0.3$ . The saddle-focus (NSF) at $\eta _{2sf} = \eta _1 \eta _3$ is indicated by an o, and the saddle node (SNB) by a +. Tipping at the NSF due to drift when $\eta _2 = 1.35 - 0.025 \; t$ is shown in dashed red. We also show the regions $S^+$ and $S^-$ , and the discontinuity set $\Sigma $ at ${\cal V} = 0$ .

The critical point $\eta _{2sf}$ , indicated by an o in Figure 1, arises at a BCB [Reference Bernardo, Budd, Champneys and Kowalczyk4]. This occurs when F and S intersect with $\Sigma $ . This critical point can be obtained from (2.1) as

$$ \begin{align*} \eta_{2sf} \equiv \eta_1\eta_3. \end{align*} $$

This (NSF) bifurcation is a nonsmooth fold [Reference Bernardo, Budd, Champneys and Kowalczyk4], in which F and S co-exist if $\eta _2> \eta _{2sf}$ and neither exist if $\eta _2 < \eta _{2sf}$ . Close to $\eta _{2sf}$ , the NSF has a V-shaped form, with the focus and saddle point having the form of $({\cal T}^{\pm },{\cal V}^{\pm })$ with ${\cal V}^{\pm }$ proportional to $\pm \mu $ , where $\mu = \eta _2 - \eta _{2sf}.$

The coalescence of a saddle S with a focus F can only occur, because this is a nonsmooth system. Such bifurcations do not arise in smooth systems where an SNB necessarily indicates collision of a stable node and a saddle. The mathematical structure near $(\eta _{2sf}, 0)$ (indicated by an o in Figure 1) is substantially different from that near the smooth SNB $(\eta _{sn}, {\cal V}_{2sn})$ (indicated by a $+$ in Figure 1). In particular, at the saddle-focus NSF, the real parts of the eigenvalues of the linearizations of either of the fixed points (saddle and focus) do not drop to zero.

2.3 Smooth and nonsmooth tipping

In general, parameters are not static in climate models of this type, but rather can oscillate seasonally and have a mean that can drift over time due to climatic variation. As an example, the fresh water hosing from the melting of the Arctic ice leads to variation in $\eta _2$ of exactly this nature. To capture the impact of such parameter variation, we consider the case where $\eta _2(t)$ is time dependent, with a mean drift of rate $\epsilon $ , and a seasonal oscillation of amplitude A and frequency $\omega $ , so that

(2.3) $$ \begin{align} \eta_2(t) = \eta_{2,0} - \epsilon t + A \sin(\omega t). \end{align} $$

Thus,

(2.4) $$ \begin{align} \begin{aligned} \dot{{\cal V}} & = \eta_1 -\eta_2(t)+\eta_3({\cal T}-{\cal V})-{\cal T}-{\cal V}|{\cal V}| ,& \dot{{\cal T}} = \eta_1-{\cal T}(1+|{\cal V}|), \end{aligned} \end{align} $$

with $\eta _2(t)$ given by (2.3). As is typical for applied settings [Reference Dijkstra8], we follow certain parameter assumptions. It is frequently assumed that the salinity’s relaxation time is much longer than that of temperature, giving $\eta _3 < 1$ , which results in an SNB in (2.4) for the branch ${\cal V}>0$ . Furthermore, we take $\eta _1 = O(1)$ so that $\eta _{2c}=\eta _3\eta _1= O(1)$ and $\eta _{2c}<\eta _{2sf}$ . That is, there is a nontrivial bi-stability range for the two stable equilibria on the branches N ( ${\cal V}>0$ ) and F ( ${\cal V}<0$ ), as shown in Figure 1. Variation of a parameter (typically $\eta _2$ in (2.2) as given by (2.3)) can lead to tipping, which in the context of this study corresponds to a solution starting at the focus F (or N) that does not stay close to F (or N) but rapidly evolves to a qualitatively different state, typically to N (or F) or to a large amplitude periodic orbit.

2.3.1 Smooth B-tipping

If $\eta _2$ is slowly increased through the SNB $\eta _{2sn}$ at a constant rate $d \eta _2/dt = \epsilon \ll 1$ , then we see canonical B-tipping with a transition from a stable temperature-dominated state (N) to a stable salinity-dominated state (F). Such tipping is described by standard theory [Reference Haberman10] and occurs at the point

(2.5) $$ \begin{align} \eta_2 = \eta_{2sn} + C_{sn} \; \epsilon^{2/3}. \end{align} $$

In this expression, $C_{sn}> 0$ , so that the tipping point is delayed till after the bifurcation (with obvious climatic implications).

2.3.2 Nonsmooth B-tipping

In contrast, if $\eta _2$ is slowly decreased through the nonsmooth fold at $\eta _{2sf}$ at a constant rate $d\eta _2/dt = -\epsilon $ , then nonsmooth B-tipping occurs, with a transition from F to N. This is illustrated in Figure 1, where we can see the lag in the tipping. B-tipping close to an NSF is different in many respects from B-tipping close to an SNB, partly because the eigenvalues of the linearization of the nonsmooth system about the fixed point do not drop to zero at a border collision bifurcation. We show presently (see also the related 1D model discussed in [Reference Budd and Kuske6]) that B-tipping occurs for the nonsmooth system, in the sense that $|V| = K \gg 1 $ , at the parameter value

$$ \begin{align*} \eta_2 = \eta_{2sf} - C_{sf} \; \epsilon \log(K/\epsilon), \quad C_{sf}>0. \end{align*} $$

Observe that the expression for the lag in the tipping point scales differently from the smooth case in (2.5).

The case of B-tipping due to seasonal forcing (possibly combined with slow drift) is more subtle and, as in the 1D case, is dominated by the existence of a cyclic-fold bifurcation which can lead to an advancement of the tipping point. In [Reference Budd and Kuske6], this was explored for an NSF in a 1D simplification of (2.2) and a pattern of tipping not dissimilar to the smooth case was observed. However, in the case of the full 2D problem (2.2), there are additional effects due to the possibility of resonance between the seasonal parameter variation and the natural frequency of the focus. This can lead to complex behaviour close to the tipping point, including chaotic dynamics. This behaviour is unique to the NSF and cannot occur at a standard SNB as there is no natural frequency in the latter case.

3 The normal form of the Stommel 2-box model close to an NSF as a forced piece-wise smooth system

3.1 Overview

Our primary interest in this paper is the nature of the tipping of the system (2.4) close to the nonsmooth saddle-focus (border collision) bifurcation. Accordingly, we track the system (2.4) near the two respective fixed points. Following an expansion around each fixed point, we find that, to leading order, the Stommel 2-box model can be simplified to a piece-wise linear normal form. This is a continuous dynamical system comprising two linear systems, separated by the surface $\Sigma \equiv \{ {\cal V} =~0 \}$ across which the solution trajectory loses smoothness. We find that the properties of the linear operators of each of the two linearized equations play a critical role in the resulting dynamics. Here, we provide a summary of the linearized system and its elementary properties.

Consider the Stommel 2-box model (2.4). At the nonsmooth fold, we have ${\cal V} = 0$ and ${\cal T} = \eta _1$ . We set ${\cal V} = V \ll 1$ and ${\cal T} = \eta _1 + T$ with $T \ll 1$ , and take

$$ \begin{align*}\eta_2(t) = \eta_1\eta_3 + \theta(t) \equiv \eta_{2sf} + \theta(t),\end{align*} $$

where from (2.3), we express the parameter drift and seasonal forcing by

$$ \begin{align*} \theta(t) \equiv \mu(t) + A\sin(\omega t) \equiv \mu_0 - \epsilon t + A \sin(\omega t), \end{align*} $$

so that a static border collision bifurcation occurs when $\mu = \epsilon = A = 0.$ Furthermore, let

(3.1) $$ \begin{align} {\mathbf x} = (T,V)^T, \quad {\mathbf e}_1 = (1,0)^T, \quad {\mathbf e}_2 = (0,1)^T. \end{align} $$

The Stommel 2-box model is smooth in the respective regions

$$ \begin{align*}S^+ = \{ x\mid {\mathbf e}_2^T {\mathbf x} = V> 0 \}, \quad S^- = \{ {\mathbf x} \mid {\mathbf e}_2^T {\mathbf x} < 0 \},\end{align*} $$

and loses smoothness across the discontinuity set

$$ \begin{align*}\Sigma = \{x\mid {\mathbf e}_2^T {\mathbf x}\} = 0.\end{align*} $$

Close to the NSF, the Stommel 2-box model then simplifies to the piece-wise linear normal form

(3.2) $$ \begin{align} \dot{\mathbf x} = L^{\pm} {\mathbf x} - \theta(t) \; {\mathbf e}_2, \end{align} $$

where the linear operators

$$ \begin{align*} L^+ = \left( \begin{array}{@{}c c@{}} -1 & -\eta_1 \\ \eta_3 - 1 & -\eta_3 \\ \end{array} \right), \quad L^- = \left( \begin{array}{@{}c c@{}} -1 & \eta_1 \\ \eta_3 -1 & -\eta_3 \\ \end{array} \right),\end{align*} $$

are respectively applied in the regions $S^+$ and $S^-$ .

Hence, the (forced) Stommel 2-box model (2.4) has the normal form

(3.3) $$ \begin{align} \dot{\mathbf x} = L^{\pm} {\mathbf x} - (\mu(t) + A \sin(\omega t)) \; {\mathbf e}_2 \equiv L^{\pm} {\mathbf x} - (\mu_0 - \epsilon t + A \sin(\omega t) )\; {\mathbf e_2}. \end{align} $$

We study the normal form (3.3) in detail for the remainder of this paper, looking at the cases of:

1. (i) drift-induced tipping close to $\mu = 0$ when $A = 0$ ;

2. (ii) seasonal tipping when $\epsilon = 0$ with $\omega $ large, small and resonant; and

3. (iii) combined tipping when $\epsilon $ and A are both nonzero.

We first derive some elementary properties of this piece-wise linear system (3.2) which we will refer to in the later Sections 4, 5 and 6.

3.2 Fixed points

It is immediate that if $\mu (t) = \mu \equiv \mu _0$ is fixed, then there are fixed points, a saddle ${\mathbf S} = \mu \; {\mathbf s}$ and a focus ${\mathbf F} = \mu \; {\mathbf f}$ given by

(3.4) $$ \begin{align} {\mathbf s} = (L^+)^{-1} {\mathbf e_2} = \frac{1}{\eta_3 + \eta_1(\eta_3 - 1)} \left( \begin{array}{@{}c@{}} \eta_1 \\ -1 \end{array} \right) , \quad {\mathbf f} = (L^-)^{-1} {\mathbf e_2} = \frac{1}{\eta_3 - \eta_1(\eta_3 - 1)} \left( \begin{array}{@{}c@{}} -\eta_1 \\ -1 \end{array} \right). \end{align} $$

Observe, as expected, that if $0 < \eta _3 < 1$ and $\eta _1 \gg \eta _3$ , then $S \in S^+$ and $F \in S^-$ if and only if $\mu> 0$ , with the BCB occurring at $\mu = 0.$ The resulting phase plane of the linearized system (3.3) either side of the NSF in the case where $A = 0$ and $\mu = \pm 1$ is given in Figure 2. In this figure, the (physical) saddle is indicated by an o and the focus (when $\mu> 1$ ) is clear from the limits of the trajectories.

Figure 2 Solutions of the linearized problem (3.3) in the $(T,V)$ phase plane when $A = 0$ . (a) $\mu = -1.$ The (unphysical) saddle is shown as an o. (b) $\mu = 1$ showing the (physical) saddle as an o and the stable focus.

We observe, trivially, that in contrast to the SNB, the linearization of the system (3.2) about the fixed points is independent of the value of $\mu $ . The corresponding linear operators $L^{\pm }$ have eigenvalues $\sigma _{1,2}^{\pm }$ and eigenvectors ${\mathbf u}^{\pm }_{1,2}$ , where

$$ \begin{align*}\sigma_1^{\pm} + \sigma_2^{\pm} = -(1 + \eta_3), \quad \sigma_1^{\pm}\sigma_2^{\pm} = \eta_3 \pm \eta_1(\eta_3 - 1),\end{align*} $$

so that if $\eta _3 \ll 1$ ,

$$ \begin{align*}\sigma_1^{\pm} = -\tfrac{1}{2} + \tfrac{1}{2} \sqrt{1 \pm 4 \eta_1}, \quad \sigma_2^{\pm} = -\tfrac{1}{2} - \tfrac{1}{2} \sqrt{1 \pm 4 \eta_1}.\end{align*} $$

Hence, for the typical range of $\eta _1,\eta _3$ , we have that $\sigma _{1,2}^+$ are real with $\sigma ^+_1> 0 > \sigma ^+_2$ , and $\sigma _{1,2}^-$ are complex with negative real part. For most of this paper, we will consider the representative values of

$$ \begin{align*}\eta_1 = 3, \quad \eta_3 = 0.3.\end{align*} $$

In this case,

$$ \begin{align*}\sigma_1^+ &= 0.84\cdots, \quad \sigma_2^+ = -2.14\cdots, \quad \sigma_1^- = -0.65\cdots + 1.41\cdots \; i, \\ &\quad\sigma_2^- = -0.65\cdots - 1.41\cdots \; i.\end{align*} $$

Accordingly, we may expect resonant behaviour in the system when the forcing frequency is close to the frequency of the focus and satisfies $\omega \approx 1.4$ .

3.3 Smooth periodic orbits in $S^-$ and grazing

Now, consider the case of $\mu $ fixed and seasonal forcing with $A> 0$ . Provided $A/\mu $ is sufficiently small, the system (3.3) has a smooth and stable periodic solution lying wholly in $S^-$ . A straightforward calculation shows that the periodic orbit is given by

$$ \begin{align*} {\mathbf P^-}(t) = \mu {\mathbf f} + A ( \omega^2 I + (L^-)^2 )^{-1} \,[ \omega \cos(\omega t) {\mathbf e_2} + \sin(\omega t) L^- {\mathbf e_2} ]. \end{align*} $$

It also has a smooth unstable periodic orbit lying wholly in $S^+$ given by

(3.5) $$ \begin{align} {\mathbf P^+}(t) = \mu {\mathbf s} + A ( \omega^2 I + (L^+)^2 )^{-1} \,[ \omega \cos(\omega t) {\mathbf e_2} + \sin(\omega t) L^+ {\mathbf e_2} ]. \end{align} $$

Both orbits play a critical role in understanding tipping close to the NSF under seasonal forcing.

If $\omega \gg 1$ , then we have to leading order that the stable periodic orbit is given by

$$ \begin{align*}{\mathbf P^-}(t) = \mu {\mathbf f} + A \frac{\cos(\omega t)}{\omega} {\mathbf e}_2.\end{align*} $$

This periodic orbit lies wholly in $S^-$ provided that

(3.6) $$ \begin{align} \frac{A}{\mu} < \omega |{\mathbf e}_2^T {\mathbf f}|. \end{align} $$

Similarly, if $\omega \ll 1$ , then to leading order, the periodic orbit is given by

(3.7) $$ \begin{align} {\mathbf P^-}(t) = \mu {\mathbf f} + A \sin(\omega t) {\mathbf f}. \end{align} $$

This orbit lies wholly in $S^-$ provided that

$$ \begin{align*}\frac{A}{\mu} < 1,\end{align*} $$

and if $A/\mu = 1$ , it touches the origin.

The stable periodic orbit ${\mathbf P}^-(t)$ in $S^-$ will have its largest amplitude when $\omega $ is close to the resonant value of $\omega = \beta $ . (Note that resonance is not important in $S^+$ .)

In general, as $\mu $ is decreased from a large value (or equivalently, as A is increased from zero), then there will be a first value $\mu _g$ (equivalently $A_g$ ) when the periodic solution ${\mathbf P}^-(t)$ grazes the line $\Sigma $ . This occurs when $\max ({\mathbf e_2}^T {\mathbf P}^-(t)) = 0$ . As $\mu $ is decreased through $\mu _g$ (A is increased), the periodic orbit will persist before vanishing at a cyclic fold at $\mu = \mu _c$ . At this point, we see cyclic B-tipping or more complex behaviour. We consider this case presently. We present in Figure 3 a plot of the (grazing) value $\mu _g$ as a function of $\omega $ in the case of $A = 1, \eta _1 = 3, \eta _3 = 0.3$ . Note that $\mu _g$ takes its largest value at the resonant frequency $\omega \approx 1.4.$ We deduce that $\mu _g$ increases close to the resonant value and, for these values of $\omega $ , we expect to see significantly advanced tipping. This is quite different from tipping in the 1D map case studied in [Reference Budd and Kuske6].

Figure 3 The grazing value $\mu _g$ , when $A=1, \eta _1 = 3, \eta _3 = 0.3$ .

4 Drift-induced B-tipping

We first consider the case of drift-only-induced tipping when $A = 0$ for which $\mu (t) = \mu _0 - \epsilon t$ and we take ${\mathbf x}(0) = \mu _0 {\mathbf f}.$ The analysis is similar to the 1D case [Reference Budd and Kuske6] and is also given in [Reference Budd, Griffth and Kuske5].

To make the calculations of tipping precise, we observe that as the normal form we are considering is a piece-wise linear system, any divergence to infinity of the solution occurs only in infinite time (unlike the case of tipping close to an SNB where (owing to the quadratic form of the nonlinearity) the divergence to infinity takes place in finite time). In the case of the piece-wise linear normal form, we (as described earlier) define tipping, in this and future sections, to have occurred when

$$ \begin{align*} {\mathbf e}_2^T {\mathbf x} = V = K \end{align*} $$

for some suitably large $K> 0$ (typically, we take $K = 10$ ). A similar definition of tipping for a 1D piece-wise linear system was used in [Reference Budd and Kuske6].

The estimate (4.1) agrees well with numerical experiments for small values of $\epsilon $ . In Figure 4, we present the tipping value of $\mu $ calculated numerically (blue) and the estimate (4.1) (red) for the case of $\eta _1 = 3, \eta _3 = 0.3$ , for which $1/\sigma _1 = 1.18$ .

The form of the tipping: (i) as a function of $\mu (t)$ and (ii) in the $(T,V)$ phase plane, is shown in Figure 5. When ${\mathbf x} \in S^-$ , and ignoring rapidly decaying exponential terms, the trajectory takes the form

$$ \begin{align*}{\mathbf x}^-(t) = \mu(t) {\mathbf f} - \epsilon (L^-)^{-1} {\mathbf f}.\end{align*} $$

This trajectory intersects the discontinuity set $\Sigma $ at the time $t_1$ when $\mu (t_1) = \epsilon \; {\mathbf e}_2^T [(L^-)^{-1} {\mathbf f}]/{\mathbf e}_2^T{\mathbf f}.$ It then enters $S^+$ , where (again ignoring rapidly decaying terms) it takes the form

$$ \begin{align*}{\mathbf x}^+(t) = \mu(t) {\mathbf s} - \epsilon (L^+)^{-1} {\mathbf s} + a \; \epsilon \; e^{\sigma_1^+ (t-t_1)} \; {\mathbf u}^+_1, \end{align*} $$

where the eigenvalue $\sigma _1^+$ is as given in Section 3, and the precise value of a is given by matching ${\mathbf x}^-$ to ${\mathbf x}^+$ on $\Sigma $ . The dominant feature of this expression is the exponentially growing term parallel to the eigenvector ${\mathbf u}^+_1$ of $L^+$ (as can be seen in Figure 5). Accordingly, we predict that B-tipping occurs, to leading order in $\epsilon $ , when

(4.1) $$ \begin{align} \mu = -\frac{\epsilon}{\sigma_1} \; \log(K/\epsilon) + {\cal O}(\epsilon). \end{align} $$

Figure 4 Drift-induced tipping comparing the calculated tipping point $\mu $ (blue) with the estimate (4.1) (red).

Figure 5 (a) Drift-induced B-tipping as a function of $\mu $ for a range of values of $\epsilon = 0.2, 0.1, 0.05, 0.01$ and with $K=10$ (from left to right respectively in green, red, black and blue). (b) Tipping (blue trajectory) in the $(T,V)$ phase plane when $\epsilon = 0.05$ . In this figure, the locus of the two points $\mu {\mathbf f} \in S^-$ and $\mu {\mathbf s} \in S^+$ are shown in red.

5 High-frequency B-tipping when $\mu $ is fixed and $\omega \gg 1$

In this section, we consider high-frequency periodic forcing with no drift, so that $\mu $ is fixed, with $\epsilon = 0$ and $\omega \gg 1$ . This case was first considered in [Reference Budd, Griffth and Kuske5] and also in [Reference Budd and Kuske6], and we extend the results obtained in these papers. The main result is then given in the following proposition.

Proposition 5.1. Let the parameter $\lambda $ be defined by

$$ \begin{align*}\lambda = \frac{A}{\omega \mu}.\end{align*} $$

Then, we have the following results.

1. (i) For $\lambda < \lambda _g = |{\mathbf e}_2^T {\mathbf f}|$ , the system (3.2) has a stable periodic solution lying entirely within $S^-$ centred on the focus $\mu {\mathbf f}$ . As $\lambda $ is decreased to $\lambda = \lambda _g$ , this orbit grazes $\Sigma $ .

2. (ii) For $\lambda _g < \lambda < \lambda _c$ , the stable periodic orbit persists and lies in both $S^-$ and $S^+$ , and for $\lambda _g < \lambda = \lambda _0 < \lambda _c$ , it has a figure-of-eight form centred on $V=0$ .

3. (iii) At $\lambda = \lambda _c$ , the stable periodic orbit ceases to exist at a cyclic fold bifurcation. For $\lambda> \lambda _c$ , the solution diverges to infinity and we see tipping behaviour.

The conclusions of Proposition 5.1 are illustrated in Figure 6. In this figure, we take $\omega = 10$ and plot the stable periodic orbit for different values of $\lambda < \lambda _c$ and an orbit (left-most) for $\lambda $ just greater than $\lambda _c$ , which is diverging to infinity. We see the transition from a near-circular orbit wholly in $S^-$ to a figure-of-eight orbit for $\lambda $ slightly smaller than $\lambda _c$ , and then to a solution which diverges to infinity (tipping behaviour) as $\lambda $ increases.

Figure 6 The behaviour under high-frequency periodic forcing with $\omega = 10$ and $\eta _1 = 3, \eta _3 = 0.3.$ Parameter values: $\lambda _g = 0.416, \lambda _0 = 0.748$ and $\lambda _c = 0.767$ . This shows the stable periodic orbits that exist for $\lambda = 0.4, 0.5,0.6,0.75,0.76$ (respectively cyan, blue, red, orange, mauve). Observe that the centre of the periodic orbit shifts to the left when it lies on both sides of $\Sigma $ and in the case of the symmetric figure-of-eight orbit, it is centred on $V = 0, T = -1/(1 - \eta _3).$ For larger $\lambda = 0.77> \lambda _c$ , the stable periodic orbit ceases to exist and, instead, the solution diverges to infinity (tipping behaviour). This behaviour is seen in the left-most orbit in the figure. The discontinuity set $V = 0$ is indicated.

We now derive the results in this proposition. Result (i) follows immediately from the results described in Section 3.2.

To derive the other results for large $\omega $ , we set ${\mathbf x} = \mu {\mathbf z}$ , $ \lambda = A/(\omega \mu )$ and $\tau = \omega t$ . Under rescaling, the linear normal form of the Stommel 2-box model becomes

$$ \begin{align*} \dot{{\mathbf z}} = L^{\pm} \; {\mathbf z}/\omega - \lambda \sin(\tau) {\mathbf e_2} - {\mathbf e_2}/\omega. \end{align*} $$

We assume that $\lambda = {\cal O}(1)$ and set

$$ \begin{align*}{\mathbf z} = {\mathbf z}_0 + {\mathbf z}_1/\omega + \cdots.\end{align*} $$

We then have, for some constant vector ${\mathbf C} = [C_1,C_2]^T$ (to be estimated at the next level of the asymptotics),

(5.1) $$ \begin{align} {\mathbf z}_0(\tau) = {\mathbf C} + \lambda \cos (\tau) {\mathbf e_2}. \end{align} $$

Hence, ${\mathbf z}_0(\tau )$ is a periodic orbit centred on ${\mathbf C}$ . Observe that in this limit, ${\mathbf z}_0(\tau )$ only oscillates in the V direction parallel to ${\mathbf e}_2$ . If $\lambda < |C_2|$ , this periodic orbit lies wholly in $S^-$ .

At the next level of asymptotics,

(5.2) $$ \begin{align} \dot{\mathbf z}_1 = L^{\pm} {\mathbf z}_0 - {\mathbf e_2}. \end{align} $$

This differential equation admits a periodic solution provided that

(5.3) $$ \begin{align} {\mathbf I} \equiv \int_0^{2\pi} L^{\pm} {\mathbf z}_0 \; d\tau = 2 \pi {\mathbf e_2}. \end{align} $$

We use this identity to determine ${\mathbf C}$ . It is immediate that ${\mathbf e_2}^T {\mathbf z}_0$ = $C_2 + \lambda \cos (\tau )$ . Hence, the choice of $\pm $ in the integrand (5.3) depends only upon $C_2$ . The following results are immediate:

(5.4) $$ \begin{align} {\mathbf e_2}^T L^{\pm} {\mathbf e_2} = -\eta_3, \quad {\mathbf e_1}^T L^{\pm} {\mathbf e_2} = \mp \eta_1, \quad L^{\pm} {\mathbf e_1} = (-1, \eta_3-1), \end{align} $$

where ${\mathbf e}_1$ and ${\mathbf e}_2$ are as defined in (3.1). Multiplying ${\mathbf I}$ by ${\mathbf e_2}^T$ and applying (5.4) gives

$$ \begin{align*} (\eta_3 - 1) C_1 - \eta_3 C_2 = 1. \end{align*} $$

We now make some estimates for $C_1$ and $C_2$ which allow us to deduce the existence of the cyclic fold at $\lambda = \lambda _c$ .

Proof. (i) To obtain this identity, we multiply ${\mathbf I}$ by ${\mathbf e_1}^T$ . This gives

$$ \begin{align*}-{2 \pi} C_1 + \int_0^{2 \pi} \mp \eta_1 (C_2 + \lambda \cos(\tau)) \, d \tau = 0.\end{align*} $$

However, by definition, the sign of $C_2 + \lambda \cos (\tau )$ is precisely the opposite of the sign of $\mp $ . Hence, the result (i). Results (ii) and (iii) follow by direct quadrature.

We now consider the effect of increasing $\lambda $ from zero to a maximum value at which there is a cyclic-fold bifurcation.

If $\lambda < |C_2|$ , the periodic orbit lies only in $S^-$ , and from the earlier results (3.4), (3.6),

$$ \begin{gather*} C_1 = -\eta_1/(\eta_3 - \eta_1 (\eta_3 - 1)), \quad C_2 = -1/(\eta_3 - \eta_1 (\eta_3 - 1)) < 0, \\ \lambda_g = 1/(\eta_3 - \eta_1 (\eta_3 - 1)). \end{gather*} $$

Now, consider the case of $\lambda> \lambda _g$ when the periodic orbit lies both in $S^+$ and $S^-$ . Combining the results above, we have that $C_2$ satisfies the identity

(5.5) $$ \begin{align} -\frac{2\eta_1(\eta_3 - 1)}{\pi} \Big[ \sqrt{\lambda^2 - C_2^2} + C_2 \sin^{-1}(C_2/\lambda) \Big] - \eta_3 C_2 = 1. \end{align} $$

The solution $C_2$ of this together with $C_1$ is plotted in Figure 7. We see that this has a (cyclic-)fold bifurcation with $\lambda $ taking its largest value $\lambda _c$ when $dC_2/d\lambda = 0$ for a small value of $C_2> 0$ . For the case plotted of $A=1, \eta _1 = 3,\eta _3=0.3$ , we have $(\lambda _c,C_2) = (0.7672,0.17).$

Figure 7 The values of $C_1$ (red) and $C_2$ (blue) as a function of $\lambda $ when $\eta _1 = 3$ and $\eta _3 = 0.3$ . The cyclic fold is clearly visible and occurs when $(\lambda _c,C_2) = (0.7672,0.17)$ .

Proof. Result (i) follows by inspection. To establish result (ii), we first assume that $C_2$ is small. We then approximate (5.5) by the quadratic equation

$$ \begin{align*}-\frac{2 \eta_1 (\eta_3 - 1)}{\pi} \bigg( \lambda + \frac{C_2^2}{2 \lambda} \bigg) - \eta_3 C_2 = 1.\end{align*} $$

Hence, in this limit,

$$ \begin{align*} \frac{1 + \eta_3 C_2}{(1 - \eta_3)} = \frac{2 \eta_1}{\pi} \bigg( \lambda + \frac{C_2^2}{2 \lambda} \bigg). \end{align*} $$

This quadratic equation has a fold bifurcation at $\lambda = \lambda _c$ when its discriminant vanishes. This gives the expressions for $\lambda _c$ and $C_2$ in (5.6). Note that if $\eta _3 \ll 1$ , then $C_2 \ll 1$ , justifying the assumptions made.

We conclude that the periodic solution $z_0$ vanishes in a cyclic-fold bifurcation for $\lambda = \lambda _c$ . For $\lambda> \lambda _c$ in this high-frequency limit, we do not have a stable periodic solution and (cyclic) B-tipping occurs. Equivalently, this tipping occurs if $A> A_c$ or if $\mu < \mu _c$ .

To complete the calculation of the form of the periodic solution, we explain the figure-of-eight shape of the orbit which occurs when $\lambda = \lambda _0$ and is shown in Figure 6. The first-order contribution to the periodic orbit is given by ${\mathbf z_1}(\tau )$ which satisfies the ordinary differential equation (ODE) (5.2) using (5.1),

$$ \begin{align*}\dot{\mathbf z}_1 = L^{\pm} (C_1 {\mathbf e_1} + (C_2 + \lambda \cos(\tau) {\mathbf e_2)}) - {\mathbf e_2}. \end{align*} $$

Let ${\mathbf z}_1 = (p,q)^T.$ By taking projections, we find that

$$ \begin{align*} \dot{p} = -C_1 -\eta_1 |C_2 + \lambda \cos(\tau)|, \quad \dot{q} = (\eta_3 - 1) C_1 -\eta_3 C_2 -1 -\eta_3 \lambda \cos(\tau). \end{align*} $$

As $(\eta _3-1) C_1 - \eta _3 C_2 - 1 = 0$ , we see that q simply oscillates in phase with ${\mathbf z}_0$ . Now, consider p. If $|C_2|$ is sufficiently large, then $C_2 - \lambda \cos (\tau )$ takes one sign only. Hence, p has the same frequency of oscillation as ${\mathbf z_0}$ . The orbits in this case will be elliptical, as we can see in the numerical calculation presented in Figure 6. In contrast, if $|C_2|$ is small and, in particular, if $\lambda = \lambda _0, C_2 = 0$ , the function $|\lambda \cos (\tau )|$ is a “rectified cosine wave”. This has period $\pi $ , half that of the function $\cos ({\tau })$ . Indeed, to a first Fourier mode approximation,

$$ \begin{align*}\lvert\cos(\tau)\rvert \approx \cos^2(\tau) = \tfrac{1}{2} (1 + \cos(2 \tau)). \end{align*} $$

Thus, in this case, p has half the period of ${\mathbf z_0}$ , which gives the figure-of-eight curve observed in this case in Figure 6.

6 Low-frequency periodic forcing with no drift

Suppose now that $\omega \ll 1$ . We have shown in Section 3.2 (3.7) that in the case of small $\omega $ and for sufficiently large $\mu $ , the orbit in $S^-$ takes the leading-order form

$$ \begin{align*}P^-(t) = {\mathbf x}(t) = [\mu + A \sin(\omega t) ] \; {\mathbf f}.\end{align*} $$

As $\mu $ decreases, at $\mu = \mu _g \approx A$ , this orbit touches $\Sigma $ close to the origin. For smaller values of $\mu $ , this trajectory enters the region $S^+$ , where it becomes subject to the exponential growth arising from the positive eigenvalue of the linear operator $L^+$ . Then, for such a periodic orbit to exist in this range, it must spend a short time in $S^+$ , limiting the effect of the exponential growth, before returning to the region $S^-$ . At a value $\mu _c < \mu _g$ , the periodic orbit loses stability at a cyclic-fold.

To pursue the construction of the periodic orbit, we consider first a scaling argument. Setting $\tau = \omega t$ with $ \omega \ll 1$ gives

$$ \begin{align*}\omega \; {\mathbf x}_{\tau} = L^{\pm} {\mathbf x} -(\mu + A \sin(\tau)) {\mathbf e}_2.\end{align*} $$

To leading order in (small) $\omega $ , noting from Section 3.1 that ${\mathbf s} = (L^+)^{-1} {\mathbf e}_2$ and ${\mathbf f} = (L^-)^{-1} {\mathbf e}_2$ , we then have

$$ \begin{align*}{\mathbf x} = \mu(1 + \sin(\tau)) {\mathbf f} \quad \mbox{or} \quad {\mathbf x} = \mu(1 + \sin(\tau)) {\mathbf s}.\end{align*} $$

Hence, the orbit will be “L-shaped” with two sections meeting close to the origin, a longer one in $S^-$ parallel to ${\mathbf f}$ and a shorter one in $S^+$ parallel to ${\mathbf s}$ . Examples of such periodic orbits, when $A=1, \eta _1 = 3, \eta _3 = 0.3$ with $\omega = 0.1, \mu = 1.008$ and $\omega = 0.25, \mu = 1.007$ , are shown in Figure 8, where the dynamics in the respective regions $S^-$ and $S^+$ is close to being parallel to the vectors ${\mathbf f}$ and ${\mathbf s}$ .

Figure 8 Periodic orbits for small $\omega $ : (a) $\omega = 0.1, \mu = 1.0008$ ; (b) $\omega = 0.25, \mu = 1.007$ . The two fixed points $\mu {\mathbf f}$ and $\mu {\mathbf s}$ are shown circled. Observe the broadening of the orbit as $\omega $ increases. The discontinuity set $V = 0$ is indicated.

Figure 9 The small $\omega $ values of the grazing point $\mu _g$ (blue), the tipping point $\mu _c$ (red) and the fitted parabolic estimate $(1 + 0.1 \omega ^2)$ (dashed magenta) when $A= 1, \eta _1 = 3, \eta _3 = 0.3$ .

In the region $S^+$ , we see linear dynamics associated with the linear operator $L^+$ . This operator has eigenvalues $\sigma _1^+> 0 > \sigma _2^+$ with corresponding eigenvectors ${\mathbf u}^+_1$ and ${\mathbf u}^+_2$ . The corresponding dynamics in $S^+$ is then given by

$$ \begin{align*}{\mathbf x}(t) = \mu {\mathbf s} + \alpha_1 e^{\sigma_1^+ t} {\mathbf u}^+_1 + \alpha_2 e^{\sigma_2^+ t} {\mathbf u}^+_2 + A(\omega^2 I + (L^+)^2)^{-1} ( \omega \cos(\omega t) {\mathbf e_2} + \sin(\omega t) L^+ {\mathbf e_2} ),\end{align*} $$

so that expanding in $\omega $ , we have the following expression for the trajectory:

(6.1) $$ \begin{align} {\mathbf x}(t) = \mu(1 + (A/\mu) \sin(\tau)) {\mathbf s} + \omega \cos(\tau) (L^+)^{-1} {\mathbf s} + \alpha_1 e^{\sigma_1^+\tau/\omega} {\mathbf u}_1^+ + \alpha_2 e^{\sigma_2^+ \tau/\omega} {\mathbf u}_2^+ + {\cal O}(\omega^2).\\[-24pt]\nonumber \end{align} $$

A similar expression for the trajectory applies in the region $S^-.$

The trajectory (6.1) enters $S^+$ close to the origin and moves towards ${\mathbf s}$ , but returns to $\Sigma $ due to the effect of the term $\alpha _1 \exp (\sigma _1^+ \tau /\omega ) {\mathbf u}_1^+.$ Because of this strong exponential growth in the region $S^+$ , the orbit can only remain in this region for a scaled time of $\tau = {\cal O}(\omega )$ . As a consequence, the deviation of the orbit in the direction towards ${\mathbf s}$ is also ${\cal O}(\omega ).$ It follows that if ${\mathbf x}$ leaves and returns to the set $\Sigma $ in this time, then $\alpha _1$ and $\alpha _2$ must also both be of ${\cal O}(\omega )$ . Hence, the length and the width of the orbit in $S^+$ are ${\cal O}(\omega )$ . Similarly, the orbit is of approximate length $2 \mu $ in $S^-$ . This behaviour can be clearly seen in Figure 8 when we compare the solutions for $\omega = 0.1$ and $\omega = 0.25$ .

In Figure 9, we present numerical experiments for small $\omega < 0.5$ in which we calculate the value of $\mu _c$ at which the periodic orbit loses stability, and compare this value with the grazing point $\mu _g$ and a fitted parabola. In this and later calculations of $\mu _c$ , we calculate the stable periodic orbit at $\mu _g$ numerically (using the Matlab routine ode15s) as the $\omega $ -limit set of the trajectory starting from the focus ${\mathbf F}$ . We then systematically reduce $\mu $ , following the orbit, until it loses stability at $\mu _c$ .

From examining this figure, we conclude that if $\omega < 0.35$ , there is a reasonable fit to a parabola so that

$$ \begin{align*} \mu_c \approx A(1 + C \; \omega^2), \quad C> 0, \end{align*} $$

where for $\eta _1 = 3, \eta _3 = 0.3$ , we have $C \approx 0.1.$ For $\omega> 0.35$ , the values of $\mu _c$ depart from the fitted parabolic curve, as the effect of resonance at $\omega \approx 1$ becomes more important.

7 Near-resonant periodic forcing with no drift

Having studied the periodic orbit for the large and small $\omega $ limits, we now look at forcing frequencies in the range of $\omega = {\cal O}(1)$ . Of special interest are values close to the natural frequency of the focus in $S^-$ , which for the canonical values of the parameters is given by $\omega = 1.4$ . The resonance linked to this leads to significantly more complex behaviour than the cases of small- and large- $\omega $ studies earlier, with a more complicated path to tipping as $\mu $ is decreased from $\mu _g$ . We see the “simple” orbit at $\mu = \mu _g$ develop additional structure through (possibly a sequence of) bifurcations, leading to (possibly) chaotic behaviour. The chaotic attractor then loses stability as $\mu $ decreases, giving rise to tipping.

To illustrate this, we fix $\omega = 1$ and consider a range of orbits as we decrease $\mu $ towards a tipping value estimated to be $\mu \equiv \mu _c= 1.205.$ These orbits are shown in Figure 10. We see here very different behaviour from before. In particular, as $\mu $ decreases, the “simple” periodic orbit does not vanish at a cyclic fold bifurcation. Instead, it has a period-doubling bifurcation at $\mu \approx 1.27$ evolving, through more bifurcations and the creation of extra loops, into a (stable) chaotic orbit as $\mu $ decreases. This chaotic orbit then ceases to exist at $\mu _c$ leading to tipping. The corresponding bifurcation diagram is shown in Figure 11(a) for $\omega = 1$ . This diagram is obtained for each $\mu $ by evolving the orbit forward till it reaches its $\omega $ -limit set and then plotting V at time intervals $t_n = 2n \pi /\omega $ . In Figure 11(b), we contrast this with a similar bifurcation diagram when $\omega = 1.2$ . In this, we also see a single period-doubling bifurcation, but the period-two orbit then becomes unstable to tipping at $\mu _c = 1.17$ without evolving into a chaotic attractor.

Figure 10 $\omega = 1$ (a–f) $\mu = 1.35,1.29,1.27,1.25,1.23,1.21$ , respectively, showing the evolution of a simple periodic orbit to a chaotic one as $\mu $ decreases.

Figure 11 Bifurcation diagram of the $omega$ -limit set of the sampled orbits as a function of $\mu $ (a) $\omega = 1$ showing the period doubling bifurcation and evolution to chaos before tipping, (b) $\omega = 1.2$ where we only see a single period-doubling bifurcation before tipping.

Because of the complexity of the behaviour close to resonance, the calculation of the location of the tipping value $\mu _c$ (so that tipping occurs if $\mu < \mu _c$ ) of $\mu $ is harder in this case. For example, we must distinguish between a (possibly) chaotic orbit tipping and tipping with a chaotic transient. Numerical plots of $\mu _c$ (obtained using the method described earlier, with $A=1,\eta _1 = 3, \eta _3 = 0.1$ ) are given in Figure 12 for the wider range of $\omega \in [0,10]$ , along with the grazing value $\mu _g$ and the large, and small, $\omega $ and estimates for $\mu _c$ . We see that $\mu _c$ agrees closely with the high-frequency estimate for values of $\omega> 3.$ As $\omega $ decreases, we can see the effect of resonance, leading to a peak in the value of $\mu _c$ at $\omega \approx 1$ . It is interesting that whilst the curve of the values of $\mu _c$ shows a resonant peak, the amplitude of this peak is less marked than the peak of the curve of the grazing values $\mu _g$ and occurs at a smaller value of $\omega $ . No such peak is seen in the case of the 1D map studied in [Reference Budd and Kuske6]. The lack of smoothness of the curve around the resonance point arises from the difficulties in exactly locating $\mu _c$ remarked on above. For values of $\mu> \mu _c$ , but close to $\mu _c$ in the “resonant region”, we observe the complex behaviour described above arising from a period-doubling bifurcation for some $\mu _{pd}$ with $\mu _c < \mu _{pd} < \mu _g.$

Figure 12 The values of the tipping point $\mu _c$ (in the case of zero drift) computed from simulations in red. These values are compared with the grazing point $\mu _g$ in blue, the high-frequency estimate $1/(\omega \lambda _c)$ for $\mu _c$ in dashed black and the low-frequency estimate $1 + 0.1 \; \omega ^2$ in dashed magenta. The impact of resonance is apparent in the sharp peak of the curve $\mu _g$ at $\omega \approx 1.4$ and the less pronounced peak of the curve $\mu _c$ close to $\omega = 1.4$ .

8 Slow drift and oscillatory forcing

We now consider the case of tipping under the combination of slow drift, so that $\mu (t) = \mu _0 - \epsilon t$ , $\epsilon> 0$ combined with oscillatory forcing of amplitude $A> 0$ . In all cases, we start the evolution from the point ${\mathbf x}(0) = \mu _0 {\mathbf f}$ and define tipping to occur when $V = {\mathbf e_2}^T {\mathbf x}(t) \equiv K = 10 $ with $\mu (t) \equiv \mu _{TP}$ at this point.

We consider two cases (respectively Case A and Case B). In Case A, we take $\mu _0 = 2$ so that at $t = 0$ (as can be seen from Figure 12), we have $\mu (0)> \mu _c$ for all values of $\omega $ . As t increases, we then see a transition from a stable solution to one which experiences tipping. In Case B, we consider the case of $\mu _0 = 1$ . In this case, (again with reference to Figure 12), we have $\mu (0) < \mu _c$ for some values of $\omega ,$ so that we immediately enter a tipping regime in this instance.

Case A: $\mu _0 = 2$ .

The form of tipping is presented in Figure 13(a,b) in the cases of $\epsilon = 0.1, 0.01$ and a variety of values of $\omega $ including small, large and values close to the resonant frequency of $\omega = 1.4$ . We can see that the value of $\mu _{TP}$ increases in the region close to the resonant values of $\omega = 1.4$ (in a similar manner to $\mu _c$ ), although this effect is combined with abrupt change in the curve (due to grazing which we describe shortly) as $\omega $ increases above one.

Figure 13 $\mu _0 = 2$ . The trajectory of $V(t)$ at tipping: (a) drift $\epsilon = 0.1 \omega = 10,5,2,1.5,1,0.25$ (blue, red, black dashed, magenta, black, green); (b) drift $\epsilon = 0.01 \omega = 10,5,2,1.5,1,0.25$ (blue, red, black dashed, magenta, black, green).

In Figure 14(a), we present the tipping value $\mu _{TP}$ as a function of $\omega $ for $\epsilon = 0.005, 0.01, 0.05, 0.1$ . Also plotted (in red) is the location of the tipping point $\mu _c$ when $\epsilon = 0$ . We note that $\mu _{TP} < \mu _c$ in all cases and that $\mu _{TP} \to \mu _c$ as $\epsilon \to 0.$ We can see two things from this picture. The first is that (as expected and seen in other related systems) increasing $\epsilon $ delays tipping. Second, we see that for the smaller values of $\epsilon = 0.005, 0.01$ , resonance has an impact on tipping. Observe that tipping is advanced for all of the values of $\omega \in [1,2]$ . Similar complexity is revealed when we plot the tipping value of $\mu $ as a function of the drift rate $\epsilon $ , see Figure 14(b).

Figure 14 $\mu _0 = 2$ . Tipping value of $\mu _{TP}$ : (a) as a function of $\omega $ with drift $\epsilon = 0.005, 0.01, 0.05, 0.1$ (green, blue, black, and magenta), showing the significant effects of resonance. Also plotted are the tipping values $\mu _c$ for $\epsilon = 0$ (red) computed earlier; (b) as a function of the drift rate $\epsilon $ when $\omega = 0.5,1,1.4,2,5$ (magenta, black, blue, dashed black, green).

For many values of $\epsilon $ and $\omega $ , particularly close to the resonant value of $\omega \approx 1.4$ , we see abrupt changes in $\mu _{TP}$ as parameters vary. This makes the tipping point $\mu _{TP}$ hard to predict. For example, if $\epsilon = 0.05, 0.1$ , then these jumps occur respectively if $\omega \equiv \omega _G \approx 1.2$ and $\omega \equiv \omega _G \approx 1.415$ . This effect is a result of the grazing phenomena identified for the 1D system in [Reference Budd and Kuske6] and [Reference Zhu, Kuske and Erneux26], which is then enhanced by resonance. As a partial explanation of this phenomenon, we note that if $\mu (t) = \mu _0 - \epsilon t$ , then there is a (unique up to terms with a decaying exponential) nonexponentially growing solution in $S^+$ given by

$$ \begin{align*}{\mathbf x_G}(t) = \mu(t) {\mathbf s} - \epsilon (L^+)^{-1} {\mathbf s} + {\mathbf P}^+(t),\end{align*} $$

where the periodic orbit ${\mathbf P}^+(t)$ is given by (3.5). Suppose that $\epsilon $ is fixed. The shift in the tipping point as $\omega $ varies through $\omega _G$ arises when the orbit starting from ${\mathbf x(0)}$ grazes the orbit ${\mathbf x}_G$ , so that if $\omega < \omega _G$ , it stays in $S^+$ and tips early, or if $\omega> \omega _G$ , it returns briefly to $S^-$ and tips significantly later. The time-dependent orbits ${\mathbf x}(t)$ for $\omega $ close to $\omega _G$ are shown in Figure 15(a) and the orbits in the phase plane in Figure 15(b). Note that the red orbit has an extra “loop” briefly re-entering $S^-$ .

Case B: $\mu _0 = 1$ .

The form of tipping is presented in Figure 16(a,b) in the cases of $\epsilon = 0.1, 0.01$ and a variety of values of $\omega $ . We can immediately see from these figures that tipping is advanced for a range of forcing frequencies close to the resonant values of $\omega = 1.4$ . This is because in this case, we, in general, have $\mu < \mu _c$ so that the trajectory in this case immediately enters a tipping state. In Figure 17(a), we present the tipping point $\mu _{TP}$ as a function of $\omega $ for $\epsilon = 0.05.01, 0.05, 0.1$ , together with the location of $\mu _c$ , and in Figure 17(b) the value of $\mu _{TP}$ as a function of $\epsilon $ in this case. We can see from the plot in panel (a) of the effect of $\omega $ on $\mu _{TP}$ that resonance has a bigger effect on tipping when compared with Case B. Observe that for all of the values of $\omega $ in the “resonant region” of $\omega \in [1,2]$ that tipping is advanced, and the range of values of $\omega $ over which the resonant affects advance is broader, and flatter, than in Case A. The curve in panel (b) is also much smoother than the related curve for $\mu _0 = 2$ . It is likely that the flattening/smoothness of these curves is a consequence of taking $\mu _0 = 1$ . In [Reference Budd and Kuske6], it was shown for the 1D start with larger $\mu _0$ , that is, farther away from tipping, that the phase of the forcing plays a more significant role, which yields more fluctuations in the tipping curve as a function of $\omega $ . It is also interesting that these figures are smoother, and more monotone, than the related figures for the 1D problem for which there is no resonance.

Figure 15 $ \mu _0 = 2$ , $\epsilon = 0.1$ . Jumps in the tipping value $\mu _{TP}$ when $\omega = 1.4, 1.41, 1.42$ (blue, red, black) varies through $\omega _G$ : (a) variation with $\mu (t)$ ; (b) phase plane plotted for the cases of $\omega = 1.4$ and $\omega = 1.41$ in blue and red. These figures show the abrupt change in the tipping point close to $\omega = 1.41$ arising at a graze of the orbit with ${\mathbf x_G}(t).$ In panel (b), the blue curve tips immediately in $S^+$ whereas the red curve has an extra loop, briefly re-enters $S^-$ and tips significantly later.

Figure 16 (a) $\mu _0 = 1$ : drift $\epsilon = 0.1 \omega = 10,5,2,1.5,1,0.25$ (blue, red, black dashed, magenta, black, green). (b) Drift $\epsilon = 0.01 \omega = 10,5,2,1.5,1,0.25,0.1$ (blue,red, black dashed, magenta, black, green).

Figure 17 $\mu _0=1$ : tipping value of $\mu _{TP}$ plotted (a) as a function of $\omega $ with drift $\epsilon = 0.005, 0.01, 0.05, 0.1$ (green, blue, black and magenta), showing the significant effects of resonance. Also plotted are the tipping values $\mu _c$ when $\epsilon = 0$ (red); (b) as a function of the drift rate $\epsilon $ when $\omega = 0.5,1,1.4,2,5$ (magenta, black, blue, dashed black, green).

9 Climate implications

The Stommel 2-box model we have studied is a significant simplification of the more complex models used to study tipping of the AMOC phenomenon, although it can be embedded as a part of much more complex models where tipping is observed [Reference Dijkstra8, Reference Ditlevsen and Ditlevsen9] and can give insights into these. In [Reference Budd and Kuske6], a summary of the climatic implications of tipping in the simplified 1D Stommel box model is given, together with related studies of the Rooth model of thermohaline circulation (THC) [Reference Lucarini and Stone14, Reference Rooth20] (which is a piecewise smooth 3-box model with a series of switches for temperature T and salinity S). Some of these are relevant to the current study which considers the tipping dynamics close to the “off state” close to the nonsmooth fold. In particular, the first conclusion that the lack of smoothness of the system at the NSF advances tipping when the parameters are subject to slow drift. Similar results are found in [Reference Wood, Rodriguez, Smith, Jackson and Hawkins25]. These results are all similar to those discussed in this paper and the statements in these papers about very large bi-stability regions in the THC models are related to features in our study. The results in this paper also extend some of these earlier observations by including the oscillatory form of the forcing and by doing so, allow the low-dimensional models to capture variability in tipping (as discussed in [Reference Ditlevsen and Ditlevsen9]). The key observation (shared with the 1D model) is that oscillatory forcing terms can advance tipping (at both an NSF and at an SNB) and that oscillatory forcing combined with drift can lead to unpredictable tipping behaviour, with large changes in the tipping time as parameters (such as the drift rate and/or the forcing frequency) are varied. This has clear implications in the predictability of tipping in climate models. A significant conclusion from the study of the tipping close to the nonsmooth focus (NSF) bifurcation of the Stommel 2-box model is that the natural frequency of the focus plays an important role in tipping. In particular, if the oscillatory forcing frequency is close to resonating with this, then we see both much more complex behaviour close to tipping (such as chaotic orbits) and also significantly advanced tipping. Whilst the analysis here is only given for the 2D model, we may well expect to have an NSF embedded as a low-dimensional reduction of a bifurcation of similar nonsmooth higher dimensional problems and with various forcings of different frequencies. This then has the potential for further resonance phenomena occurring due to the existence of the focus with again complex dynamics and advanced tipping with clear climate implications. We repeat that by concentrating on the dynamics close to the NSF where we tip from the “off state” to the “on state”, we have only studied a small part of the dynamics of the Stommel 2-box or more general climate models, and that the more physically relevant (to the AMOC) tipping from the “on state” is best described by the analysis for an SNB. We emphasize that comparable resonance phenomena to those seen in this paper do not appear close to tipping at a smooth SNB, and this is a novel feature of the NSF and any physical system where an NSF might occur. However, it is important to note that transitions both at the SNB and at the NSF contribute to the hysteresis window of stable behaviour seen when moving between states, for example, effecting the return to the on state in a hysteresis loop, and this impacts on the full dynamics of the climate model.