Abrupt climate change as a rate-dependent cascading tipping point

We propose a conceptual model comprising a cascade of tipping points as a mechanism for past abrupt climate changes. In the model, changes in a control parameter, which could for instance be related to changes in the atmospheric circulation, induce sequential tipping of sea ice cover and the ocean's meridional overturning circulation. The ocean component, represented by the well-known Stommel box model, is shown to display so-called rate-induced tipping. Here, an abrupt resurgence of the overturning circulation is induced before a bifurcation point is reached due to the fast rate of change of the sea ice. Because of the multi-scale nature of the climate system, this type of tipping cascade may also be a risk concerning future global warming. The relatively short timescales involved make it challenging to detect these tipping points from observations. However, with our conceptual model we find that there can be a significant delay in the tipping because the system is attracted by the stable manifold of a saddle during the rate-induced transition before escaping towards the undesired state. This opens up the possibility for an early warning of the impending abrupt transition via detection of the changing linear stability in the vicinity of the saddle. To do so, we propose estimating the Jacobian from the noisy time series. This is shown to be a useful generic precursor to detect rate-induced tipping.

1 Introduction

Multiple elements in the Earth system are believed to be at risk of undergoing abrupt and irreversible changes in response to rising atmospheric greenhouse gas concentrations. Among others, the Arctic sea ice, the Greenland and West Antarctic ice sheets, the Amazon rainforest, and the Atlantic Meridional Overturning Circulation (AMOC) have been identified as potentially crossing such tipping points at varying levels of global warming (Lenton et al., 2008). While an abrupt decline of the Arctic sea ice is already well underway (IPCC, 2019), for a system like the AMOC it is much more uncertain if and when a tipping point will be reached. Nevertheless, it constitutes a risk that deserves attention as it has been observed across the hierarchy of climate models (Weijer et al., 2019), and there is evidence that it has occurred repeatedly during the last glacial period (Henry et al., 2016). Such past changes in the AMOC likely led to abrupt climate changes known as Dansgaard–Oeschger (DO) events (Dansgaard et al., 1993). These are the most significant instances of large-scale climate change in the past, but the underlying mechanisms remain debated.

Mathematically, tipping points are typically understood as a transition from one stable attractor of the system to another. Most often, this transition is associated with a bifurcation or attractor crisis, where a system state loses stability as a critical threshold in a control parameter is crossed, leading to tipping to another attractor (bifurcation tipping). However, tipping can occur also before a critical threshold is crossed. Stochastic perturbations may induce a transition to an alternative attractor (noise-induced tipping). Furthermore, some systems can fail to track their moving equilibrium and tip to another attractor while no bifurcation was crossed given there is a change in a control parameter at a high enough rate (rate-induced tipping) (Wieczorek et al., 2011; Ashwin et al., 2012).

Such rate-induced transitions can be expected to play a role in systems that are comprised of coupled components with a timescale separation. Here, changes in one component alter the conditions of another and act as a rapidly changing control parameter that could cause a rate-induced transition. This might occur in the real climate system, where a vast range of timescales is present in the atmosphere, ocean and cryosphere, and where important climate parameters, such as polar ice melt, currently display accelerating rates of change (Trusel et al., 2018; Bevis et al., 2019; The IMBIE Team, 2020). Indeed, a rate-induced collapse of the AMOC has been shown recently in a global ocean model (Lohmann and Ditlevsen, 2021). Rate-induced transitions in coupled systems are an even higher risk if one of the subsystems experiences abrupt change due to tipping. This constitutes a cascade of subsequent tipping points. Tipping cascades in coupled ecological or climate models have been considered before (Cai et al., 2016; Dekker et al., 2018; Rocha et al., 2018; Klose et al., 2020; Wunderling et al., 2020). However, cascades where subsystems permit rate-induced tipping have not been studied yet.

Here we explore such a scenario with a conceptual sea ice–ocean model. The model describes the influence of changing polar sea ice cover on the AMOC and features the possibility of a rate-induced resurgence of the AMOC. While an AMOC resurgence is not an issue for contemporary climate change, it plays an important role in past abrupt climate changes and DO events in particular, where it is thought to be responsible for the transitions from cold (so-called stadial) periods to prolonged episodes of milder (interstadial) conditions during the last glacial period. It is still unknown what drove these transitions and the associated resurgences of the AMOC. In climate models, an abrupt collapse of the AMOC can be induced by sudden discharges of freshwater into the North Atlantic, which is a phenomenon known to have occurred in the past (Heinrich, 1988). Similar events of sudden “removal” of freshwater that potentially lead to an abrupt resurgence of the AMOC are less well known. Instead, we consider changes in atmosphere–ocean heat exchange as driver of the AMOC resurgence. These could result from abrupt changes in sea ice cover, which in turn could be driven by changing atmospheric wind stress. The potential of rapid sea ice changes to advance the abrupt DO warming events is well established (Li et al., 2005; Dokken et al., 2013; Vettoretti and Peltier, 2016; Sadatzki et al., 2019) and has been translated into a number of conceptual models before. Gottwald proposes a model with sea ice as an intermittent thermal insulator to the polar ocean, forced by a chaotic (quasi-stochastic) atmospheric component, extremes of which can trigger temporary excursions of the ocean circulation (Gottwald, 2021). While we include a stochastic forcing, the main cause of the abrupt transitions in our model is a deterministic underlying parameter shift. A different conceptual model by Boers et al. (2018) considers sea ice and an ice shelf coupled to an ocean box model, where the sea ice evolves due to a prescribed piecewise-linear feedback, leading to self-sustained oscillations. The mechanism proposed here is different in that it involves a cascade: a tipping of the sea ice cover due to slowly changing climatic conditions leads to a rate-induced tipping of the ocean circulation as a consequence of the rapid increase in ocean heat loss.

Several lines of evidence from proxy data and climate model simulations motivate such a sequence of events. Zhang et al. (2014) showed model simulations with abrupt climate changes similar to DO events by gradually varying the Northern Hemisphere ice sheet topography, which led to shifts in the atmospheric circulation that altered the wind-driven export of sea ice. This eventually led to an abrupt decrease in North Atlantic sea ice cover and a resurgence of the AMOC. Kleppin and co-workers reported spontaneous transitions of the AMOC that were triggered by the stochastic atmospheric forcing and subsequent changes in North Atlantic sea ice (Kleppin et al., 2015). Ice core data indicate that abrupt shifts in the sea ice extent at the onset of DO events were preceded by shifts in atmospheric circulation by about a decade (Erhardt et al., 2019). Furthermore, there is evidence for gradual trends leading up to the abrupt shifts in both sea ice and atmospheric circulation proxies, indicating an underlying parameter shift that might be mutually expressed in sea ice and atmosphere (Lohmann, 2019; Sadatzki et al., 2019).

Besides illustrating a new mechanism for abrupt climate change, the conceptual model proposed here gives some interesting insight into dynamical phenomena in systems combining time-dependent and stochastic forcing. We find that the ocean component of our model (the well-known Stommel box model) displays rate-induced tipping in what could be called a “soft” tipping point. Here, due to a non-smooth fold bifurcation, tipping always occurs before the bifurcation point is reached, even if the rate of change in the parameter shift is arbitrarily slow. Further, the rate-induced transition involves attraction by the stable manifold of a saddle, which can lead to a significant delay of the tipping under stochastic forcing. Based on this, we propose an early warning indicator to detect rate-induced tipping; so far only early warning signals specific to bifurcation tipping are known (Held and Kleinen, 2004; Dakos et al., 2008; Scheffer et al., 2009, 2012).

The paper is structured as follows. In Sect. 2 the coupled conceptual model is presented. We then show rate-induced tipping of the ocean component (the Stommel box model) in the deterministic and stochastic case in Sects. 3.1 and 3.2, respectively. Thereafter, the cascading dynamics of the coupled model are presented (Sect. 3.3). Early-warning signals for the cascade, as well as for the rate-induced tipping, are investigated in Sects. 3.4 and 3.5. The results are discussed in Sect. 4, and our conclusions are given in Sect. 5.

Figure 1 (a) Bifurcation diagram of the Stommel box model with η₁ as control parameter, η₂=1.0 and η₃=0.3. Solid lines indicate branches of stable fixed points, whereas dotted lines indicate unstable fixed (or saddle) points. (b) Bifurcation diagram with η₂ as control parameter, η₁=3.0 and η₃=0.3. (c) Dependence of bi-stability on η₃ (η₁=3.0). The individual bifurcation diagrams with η₂ as control parameter are shown with decreasing bistability interval as η₃ is increased from 0.1 up to 0.75.

2 Model

2.1 Ocean component: Stommel's 1961 box model

We consider the Stommel box model of the Atlantic thermohaline circulation (Stommel, 1961) with added noise to represent variations in the atmospheric forcing on short timescales. The model assumes the overturning flow ψ in between well-mixed polar and equatorial ocean basins to be proportional to the density difference

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }\mathit{\psi }\propto ({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{e}})=\left[{\mathit{\alpha }}_T\right(T_{\mathrm{e}}-T_{\mathrm{p}})-{\mathit{\alpha }}_S(S_{\mathrm{e}}-S_{\mathrm{p}}\left)\right],\end{array}$$

where the density is given by the equation of state of seawater

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }{\mathit{\rho }}_{\text{e,p}}={\mathit{\rho }}_{\mathrm{0}}\left[\mathrm{1}-{\mathit{\alpha }}_T(T_{\text{e,p}}-T_{\mathrm{0}})+{\mathit{\alpha }}_S(S_{\text{e,p}}-S_{\mathrm{0}})\right],\end{array}$$

with some reference densities, temperatures and salinities represented by ρ₀, T₀ and S₀, respectively. The two model variables represent the dimensionless temperature difference $T={\mathit{\alpha }}_T(T_{\mathrm{e}}-T_{\mathrm{p}})$ and salinity difference $S={\mathit{\alpha }}_S(S_{\mathrm{e}}-S_{\mathrm{p}})$ in between the boxes. This defines the dimensionless overturning circulation strength

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }q=T-S.\end{array}$$

Temperature and salinity in the boxes relax towards an atmospheric temperature and salinity forcing $T_{\text{e,p}}^{\mathrm{a}}$ and $S_{\text{e,p}}^{\mathrm{a}}$. The meridional difference of the forcing drives the circulation and is represented by the two parameters ${\mathit{\eta }}_{\mathrm{1}}\propto (T_{\mathrm{e}}^{\mathrm{a}}-T_{\mathrm{p}}^{\mathrm{a}})$ and ${\mathit{\eta }}_{\mathrm{2}}\propto (S_{\mathrm{e}}^{\mathrm{a}}-S_{\mathrm{p}}^{\mathrm{a}})$. A third parameter represents the timescale ratio of the temperature and salinity relaxation ${\mathit{\eta }}_{\mathrm{3}}=\frac{{\mathit{\tau }}_T}{{\mathit{\tau }}_S}$. The model is then defined by the stochastic differential equations

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }\mathrm{d}{T}_t=\left({\mathit{\eta }}_{\mathrm{1}}-T-|T-S|T\right)\mathrm{d}t+{\mathit{\sigma }}_T\mathrm{d}{W}_{T,t}\text{(5)}& {\displaystyle }& {\displaystyle }\mathrm{d}{S}_t=\left({\mathit{\eta }}_{\mathrm{2}}-{\mathit{\eta }}_{\mathrm{3}}S-|T-S|S\right)\mathrm{d}t+{\mathit{\sigma }}_S\mathrm{d}{W}_{S,t},\end{array}$$

with the Wiener processes W_S,t and W_T,t. Time is scaled with respect to the ocean timescale τ_T=200 years. For a more detailed derivation of the model, see Dijkstra (2008). The deterministic system features a parameter regime with two stable equilibria, which are referred to as the circulation “on” and “off” states. For the on state we have T>S, where the temperature forcing gradient dominates the opposing salinity forcing gradient and drives the circulation. The off state (S>T) corresponds to a reversed circulation, which is weaker and dominated by the salinity forcing gradient. In Fig. 1a and b we show deterministic bifurcation diagrams of q with respect to the parameters η₁ and η₂. In both cases, the on state loses stability in a regular saddle-node bifurcation, whereas the off state destabilizes in a non-smooth saddle-node bifurcation. The latter is also known as a non-smooth fold (di Bernardo et al., 2008) and is due to the fact that the Stommel model is a non-smooth dynamical system owing to the absolute value in its equations (see Sect. S3 and Fig. S3 for more detail). The existence and extent of bi-stability depends on the parameter η₃. A large timescale separation (slower salinity damping) leads to a large region of bi-stability, whereas as the salinity damping approaches the timescale of temperature damping, the bistability disappears (Fig. 1c). This is because a faster salinity damping disables the positive salt advection feedback, which gives rise to the bi-stability.

2.2 Coupled sea ice–ocean model

The ocean model is coupled to a sea ice component in the polar ocean box, which is the energy balance model described in Eisenman and Wettlaufer (2009) and Eisenman (2012), modified by neglecting the seasonal cycle and effects of the sea ice thickness. The changing sea ice cover acts to insulate the polar ocean to varying degrees from the cold atmospheric temperature forcing $T_{\mathrm{p}}^{\mathrm{a}}$, thus modulating the temperature forcing gradient ${\mathit{\eta }}_{\mathrm{1}}\propto (T_{\mathrm{e}}^{\mathrm{a}}-T_{\mathrm{p}}^{\mathrm{a}})$. A schematic of the coupled model, including model variables and important parameters, is given in Fig. 2. The deterministic sea ice component is defined (Eisenman and Wettlaufer, 2009) by

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}}=\mathrm{\Delta }\tanh \left({\displaystyle \frac{I}{h}}\right)+\left[R_{\mathrm{0}}\mathrm{\Theta }\left(I\right)-B\right]I+L-F-\mathrm{1}+R,\end{array}$$

with the Heaviside step function Θ(⋅). Time is scaled with respect to τ_I=1 year. The parameters and their values are described in Table 1. While I>0 corresponds to a positive sea ice cover, I<0 represents zero sea ice cover, and the variable is instead a measure of the enthalpy of the surface ocean (Eisenman, 2012). The control parameter R models influences on the sea ice concentration due to external factors, such as export or import of sea ice into the North Atlantic via changes in wind stress. While in the climate system R is driven by slower dynamic processes, such as changes in ice sheet topography, we treat it as a control parameter. We use parameter values from Eisenman (2012), which yield a sea ice component that is bi-stable with respect to R. As seen in Fig. 3, for a range of R there exists a stable state with a positive sea ice cover (red), as well as a state with zero sea ice cover I<0 (black). This range is bounded by two saddle-node bifurcations. The stable state with sea ice cover collapses at $R=-\mathrm{0.282}$. We define the state at R=0 as the stadial state, yielding a fixed point with positive sea ice cover $I_{\mathrm{0}}^{+}\approx \mathrm{1.156}$. A slight deviation from the parameters of Eisenman (2012) is our larger value of h, which gives a more gradual albedo transition from an ice-free to an ice-covered state. Accordingly, the bifurcation diagram is more “S-shaped” instead of “Z-shaped” (see Sect. S1 and Fig. S1 in the Supplement for more details).

Figure 2 Schematic of the coupled sea ice–ocean model including model parameters and variables (bold). A description of the parameters is given in Table 1. The well-mixed polar and equatorial ocean boxes are coupled by a surface flow q, along with an identical return flow at the bottom. The ocean component is reduced to the two variables $T\propto T_{\mathrm{e}}-T_{\mathrm{p}}$ and $S\propto S_{\mathrm{e}}-S_{\mathrm{p}}$. In the polar ocean box, the sea ice cover I insulates the ocean from the cold atmospheric temperature $T_{\mathrm{p}}^{\mathrm{a}}$.

Figure 3 Bifurcation diagram of the sea ice component for parameter values as in Table 1. The solid (dotted) lines indicate stable (unstable) fixed points.

Table 1 Description of model parameters.

To model transitions from stadial to interstadial conditions in the coupled sea ice–ocean model, we consider the following mechanism. The glacial polar ocean is insulated by a high sea ice concentration from the atmospheric temperature forcing, preventing it from losing heat efficiently. As the sea ice concentration decreases, the polar ocean becomes more and more exposed to the cold atmosphere and loses heat. Thus, the sea ice variable modulates the parameter η₁, which we now define as ${\mathit{\eta }}_{\mathrm{1}}\left(I\right)={\mathit{\eta }}_{\mathrm{1}}^{\mathrm{0}}-\mathit{\kappa }\mathrm{\Theta }\left(I\right)I$, with the Heaviside function Θ(⋅) since I<0 corresponds to zero sea ice cover. Adding noise as a model of fast atmospheric perturbations (Wiener process W_I,t), this yields the following coupled equations:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{d}{I}_t=\left(\mathrm{\Delta }\tanh \left({\displaystyle \frac{I}{h}}\right)+\left[R_{\mathrm{0}}\mathrm{\Theta }\right(I)-B]I+L-F-\mathrm{1}+R\right)\mathrm{d}t\\ \text{(7)}& {\displaystyle }& {\displaystyle }+{\mathit{\sigma }}_I\mathrm{d}{W}_{I,t},\text{(8)}& {\displaystyle }& {\displaystyle \frac{{\mathit{\tau }}_T}{{\mathit{\tau }}_I}}\mathrm{d}{T}_t=\left({\mathit{\eta }}_{\mathrm{1}}^{\mathrm{0}}-\mathit{\kappa }\mathrm{\Theta }\left(I\right)\cdot I-T-|T-S|T\right)\mathrm{d}t+{\mathit{\sigma }}_T\mathrm{d}{W}_{T,t},\text{(9)}& {\displaystyle }& {\displaystyle \frac{{\mathit{\tau }}_T}{{\mathit{\tau }}_I}}\mathrm{d}{S}_t=\left({\mathit{\eta }}_{\mathrm{2}}-{\mathit{\eta }}_{\mathrm{3}}S-|T-S|S\right)\mathrm{d}t+{\mathit{\sigma }}_S\mathrm{d}{W}_{S,t}.\end{array}$$

The value of κ reflects the change in atmospheric temperature forcing when removing the sea ice cover. In this conceptual framework it can only be chosen heuristically. We can for instance assume ${\mathit{\eta }}_{\mathrm{1}}^{\mathrm{0}}=\mathrm{3.0}$ for an open ocean, and atmospheric temperature forcings in a glacial climate of 20 and −10 ^∘C in the equatorial and polar box, respectively. Full sea ice cover would limit the polar temperature forcing to 0 ^∘C, corresponding to η₁=2.0. Even if the glacial polar atmosphere were above 0 ^∘C, given that it was colder than the surface ocean, extensive sea ice cover would severely reduce heat loss to the atmosphere and thus effectively reduce η₁. Here we choose a scenario where during the stadial the sea ice reduces the atmospheric forcing from η₁=3.0 to η₁=2.65. κ is then chosen such that η₁=2.65 at the stadial fixed point $I_{\mathrm{0}}^{+}$ and η₁=3.0 for I<0, yielding $\mathit{\kappa }=\mathrm{0.35}/I_{\mathrm{0}}^{+}$. As a result, the ocean component is in the bi-stable regime for both full and zero sea ice cover. A transition from stadial to interstadial will then be captured by decreasing R from zero beyond the bifurcation point, which tips the sea ice component towards a state of I<0, while the ocean remains in the bi-stable regime.

Due to the unidirectional, linear coupling of the model, and our focus on a specific dynamical regime, we restrict our presentation of the coupled dynamics to the individual bifurcation diagrams of the sea ice component with R as control parameter and of the ocean component with η₁(I) as effective control parameter. The full bifurcation structure of the coupled model with R as the only control parameter is presented in Sect. S2 and Fig. S2.

3 Results

3.1 Rate-induced tipping and soft tipping points in the Stommel model

Here we investigate the tipping dynamics in the ocean component in the deterministic limit (${\mathit{\sigma }}_T={\mathit{\sigma }}_S=\mathrm{0}$). As noted above, there is a non-smooth fold in the Stommel model as the off state loses stability, which leads to a resurgence of the AMOC. In the bifurcation diagrams of Fig. 4 it can be seen that, both in terms of T and S, the stable fixed point (red line) moves in the same direction as the saddle point when the non-smooth fold bifurcation is approached. Thus, in a sufficiently fast parameter shift towards the fold, the saddle point can outpace the system state, which is trying to follow the moving equilibrium. This is illustrated in Fig. 4, where instantaneous parameter shifts and the corresponding movements of the system state vector in the bifurcation diagrams are indicated. When the saddle point moves past the system state, the system will tip towards the alternative stable state, which is the on circulation in our case. Thus, tipping can occur even before the bifurcations points are reached, which is known as rate-induced tipping. While in the Stommel model this can happen for both η₁ and η₂ as control parameter, it occurs for a larger range of amplitudes and rates of the parameter shift when changing η₁.

Figure 4 (a, b) Bifurcation diagram of the Stommel model equilibria in terms of the variables (a) T and (b) S as a function of η₁ as control parameter with η₂=1.0 and η₃=0.3. (c, d) The same as (a, b) but with η₂ as control parameter with η₁=3.0 and η₃=0.3. Solid lines indicate stable fixed points, whereas dotted lines indicate saddle points. The horizontal arrows indicate the movement of the system state as the control parameter is changed instantaneously within the bi-stable regime. In (a) we illustrate how the system state may track the moving equilibrium for a slow parameter shift (purple trajectory) or tip to the undesired equilibrium in a fast parameter change (blue trajectory).

To be more rigorous one has to consider the movement of the basin boundary as the control parameter is changed. The basin boundary is the stable manifold of the saddle, and it separates the basins of attractions of the on and off states, i.e., the sets of initial conditions that converge to the respective attractors. In Fig. 5 we illustrate the movement of the fixed points and basin boundary as η₁ is changed from 2.65 to 3.0. This corresponds to the scenario of a transition from stadial to interstadial sea ice cover in the coupled model, as described in Sect. 2.2. Figure 5b shows that the off fixed point before the parameter shift (open circle) lies inside the basin of attraction of the on fixed point after parameter shift (blue area). This is a sufficient condition for rate-induced tipping, which has been called basin instability (O'Keeffe and Wieczorek, 2020), since for an instantaneous parameter shift, the system would tip to the other attractor. Similarly, as the system tries to follow the moving fixed point during a sufficiently fast parameter shift, it will fail to reach the off basin (orange area) at the end of the parameter shift and tip to the on fixed point. This happens for the blue trajectory, where the parameter is ramped up linearly within 300 years. In contrast, the purple trajectory shows that tipping does not occur for a ramping duration of 500 years. For this given amplitude of the parameter shift, there is a critical rate of parameter change in between these two values.

Figure 5 Phase portraits of the Stommel model with basins of attraction and fixed points. Squares and dots indicate stable fixed points, and triangles denote saddle points. (a) Phase portrait for η₁=2.65 with several flow lines to indicate the dynamics around the saddle. The basin of attraction of the off (on) state is shaded in orange (blue). (b) Phase portrait for η₁=3.0. Two trajectories where η₁ is ramped linearly from η₁=2.65 to η₁=3.0 within 300 and 500 years are shown in blue and purple, respectively. The initial conditions $T,S=(\mathrm{2.4},\mathrm{2.5})$ are indicated by the yellow cross. Open symbols indicate the positions of the fixed points at η₁=2.65, and the black curve indicates the corresponding basin boundary from (a).

Figure 6 shows time series of q for simulations with different ramping durations. The realizations in Fig. 6a and b tip to the on attractor, while the realizations in Fig. 6c and d track the moving off equilibrium. The critical ramping duration is in between the 388.5 and 390 years employed in Fig. 6b and c. Comparing Fig. 6a to b, one observes a delay in the tipping in Fig. 6b of thousands of years. This occurs because, for close-to-critical rates, the system state passes by very closely to the saddle point, where it remains for a long time as the dynamics slow down before being ejected. The close approach of the saddle happens because the system state is attracted by the saddle's stable manifold, which is also the basin boundary. If one were to use the exact critical ramping duration, the system state would evolve precisely towards the saddle and remain there. Such trajectories are called maximum canards (O'Keeffe and Wieczorek, 2020). This behavior is also seen for trajectories that eventually track the moving equilibrium, as in Fig. 6c. It is worth noting that the attraction by the stable manifold of the saddle continues after the parameter shift is already over, as shown in the inset in Fig. 6c.

Figure 6 Time series of $q=T-S$ in the Stommel model when ramping the parameter from η₁=2.65 to η₁=3.0 at different rates. The realizations are initialized at $T,S=(\mathrm{2.4},\mathrm{2.5})$, which is close to the off fixed point at η₁=2.65. The duration of the ramping is indicated by the gray shading. The realizations in (a) and (b) with ramping durations of 300 and 388.5 years, respectively, tip from the off to the on attractor. The realizations in (c) and (d) with ramping durations of 389 and 500 years, respectively, track the moving off attractor. The on, off and saddle fixed points at η₁=3.0 are shown as horizontal lines.

Figure 7 (a) Critical ramping duration below which there is a rate-induced tipping in the Stommel model when shifting the parameter from η₁=2.65 to ${\mathit{\eta }}_{\mathrm{1}}=\mathrm{2.65}+\mathrm{\Delta }{\mathit{\eta }}_{\mathrm{1}}$. The off attractor loses stability in the bifurcation at ${\mathit{\eta }}_{\mathrm{1}}^{\text{off}}=\mathrm{3.333}$, as indicated by the dashed red line. (b) Normalized shortest distance to the basin boundary Δℬ as a function of the normalized distance to the bifurcation $\mathrm{\Delta }\mathit{\lambda }=({\mathit{\eta }}_{\mathrm{1}}^{\text{off}}-{\mathit{\eta }}_{\mathrm{1}})\cdot ({\mathit{\eta }}_{\mathrm{1}}^{\text{off}}-{\mathit{\eta }}_{\mathrm{1}}^{\text{on}}{)}^{-\mathrm{1}}$. ${\mathit{\eta }}_{\mathrm{1}}^{\text{on}}$ is the parameter value at the other saddle node bifurcation of the on state. The solid black curve shows the results of the Stommel model, and a proposed super-linear scaling is shown by the dashed curve. Also shown are results for the smooth bifurcation in the sea ice component (solid blue) and the corresponding square-root scaling (dotted).

Figure 8 (a) Probability of a rate-induced tipping in the Stommel model from the off to the on state as a function of the linear parameter ramping duration from η₁=2.65 to η₁=3.00. Different noise levels ${\mathit{\sigma }}_T={\mathit{\sigma }}_S=\mathit{\sigma }$ are considered: σ=0.01 (lightest gray curve), σ=0.02, σ=0.04, σ=0.06, σ=0.1, σ=0.2, and σ=0.4 (darkest gray curve). The dashed red line is the critical ramping duration in the deterministic system. (b) Probability distributions of the time of tipping, defined by the first crossing of q>0.1, for different noise levels. The ramping is started in year 1000 and the duration is fixed at 300 years. The dashed red line is the time of tipping in the deterministic system.

The critical ramping duration depends on the amplitude Δη₁ of the parameter shift (Fig. 7a). Rate-induced tipping becomes possible at a certain Δη₁, where the basin instability condition is first satisfied. Increasing Δη₁ then leads to a very rapid increase of the critical ramping duration D_c. Thereafter, D_c keeps increasing and actually diverges as the bifurcation is approached. This is due to the non-smoothness of the fold bifurcation, where the attractor and saddle meet in a cusp (see Sect. S3 and Fig. S3 for more detail). As a result, the attractor gets close to the basin boundary very quickly as the bifurcation is approached. This leads to a super-linear scaling of the shortest distance to the basin boundary. In Fig. 7b, we compare this to the square-root scaling of the smooth fold bifurcation in the sea ice component. In the non-smooth case, as the bifurcation is approached the basin boundary gets arbitrarily close to the attractor. Then, even very small and slow parameter increases lead to tipping. Thus, the non-smooth fold leads to what could be called a soft tipping point. In practice, there is no hard critical threshold of the parameter, but for any parameter shift at finite rate, the tipping will occur at some point prior to the bifurcation. The precise location of the tipping point depends on the trajectory of the parameter shift.

3.2 Noisy rate-induced tipping

We now consider additive noise in the ocean component, which models variations in atmospheric forcing on short timescales. In addition to the soft tipping just described, the stochastic perturbations further blur the critical threshold leading to tipping. For a given amplitude of the parameter shift, there is no longer a critical rate but a range of rates where the probability of tipping goes from 0 to 1. Figure 8a shows how this range of rates expands for increasing noise level. Note that since the system features unbounded noise, we consider finite time-tipping probabilities during a simulation time of 5000 years. Eventually, there will always occur a noise-induced transition to the on attractor, especially from the off attractor at η₁=3.0 for higher noise levels.

By introducing noise, tipping becomes a mixture of rate-induced and noise-induced transitions, since the unbounded noise allows the system to cross the basin boundaries of the deterministic system in any circumstances. Still, for low noise levels the behavior strongly resembles the deterministic case. As discussed earlier, for a ramping speed relatively close to the critical rate, the tipping involves an escape from the saddle. This behavior is robust for low noise levels, where the stochastic fluctuations cannot overcome the attraction of the stable manifold of the saddle. Thus, the system approaches the saddle before being ejected from its vicinity.

Figure 9 (a–c) Three realizations in phase space of the Stommel model with ${\mathit{\sigma }}_T={\mathit{\sigma }}_S=\mathrm{0.2}$, where η₁ is ramped from η₁=2.65 to η₁=3.00 over 300 years. The filled dot (triangle) marks the off fixed point (saddle) at η₁=3.0. The colored areas are the quasi-stationary basins of attraction at the time when their boundary is first crossed. The colored basins of attractions are given at the time of first basin crossing of the trajectories, which change color from purple to yellow. The initial conditions $T,S=(\mathrm{2.4},\mathrm{2.5})$ are indicated by the yellow cross. The locations of the saddle point (triangle) and the on fixed point at this time are shown with open symbols. The threshold $q=T-S=\mathrm{0.1}$ used to define the time of tipping is shown as a dotted line.

As the noise level is increased, there are noise-induced early tippings as well as significantly delayed tippings. In order to quantify when a tipping is “early” or “late”, we need to define the moment when the system actually tips. For the deterministic system, a sensible choice would be the time when the moving, quasi-stationary basin boundary is crossed, since this is the first moment that the system would tip in case the parameter shift would be stopped suddenly. However, for the noisy system this does not guarantee tipping, since the system may cross back to the other basin at any time. As a heuristic definition of tipping, we can instead detect the departure from the vicinity of the saddle in terms of the overturning q, as the tipping is associated with a monotonic increase of q (see Fig. 6). Thus, as tipping we define the first crossing of q=0.1, which is a slightly larger value than at the saddle to allow for some fluctuations around it. In phase space this defines a straight line.

Figure 9 shows the crossing of this threshold, as well as the basins at the time when the basin boundary is first crossed, for three different realizations with a ramping duration of 300 years and ${\mathit{\sigma }}_T={\mathit{\sigma }}_S=\mathrm{0.2}$. The time of tipping varies significantly and depends primarily on the proximity of the approach to the saddle and the subsequent time spent in its vicinity. Whereas Fig. 9b shows a realization with tipping close to the deterministic scenario, the realization in Fig. 9a leaves the stable manifold early and does not approach the saddle closely. The realization in Fig. 9c approaches the saddle very closely and remains there for a long period of time.

The tipping time distribution and its dependence on the noise level is shown in Fig. 8b. In our case of a ramping duration slightly below the critical value of the deterministic system, there are three regimes of noise levels. For low noise (σ=0.01, σ=0.02, σ=0.04 and σ=0.06 in Fig. 8b) the trajectories are very similar to the deterministic case, and it is very unlikely that the noise pushes the system closer to the saddle. Thus, the tipping time is distributed closely around the deterministic value. For intermediate noise (σ=0.1 and σ=0.2 in Fig. 8b), some early noise-assisted tippings are possible, as seen by the shift of the distribution mode towards earlier times. For other realizations there is a good chance that the noise pushes the system closer to the saddle, where it can stay for a long time (thousands of years) as the dynamics slow down before escaping. This leads to the development of a long tail in the tipping time distribution. For larger noise (σ=0.4 in Fig. 8b), even earlier noise-assisted tippings are seen, as well as some delayed tippings. However, the latter do not occur as frequently as for intermediate noise since the average residence time at the saddle is also shortened.

3.3 Cascading dynamics

We now consider the coupled model and investigate how a stadial–interstadial transition can arise as a cascading tipping of the two components. The cascade is initiated by a change in the control parameter R leading to a decrease and eventual tipping of the sea ice to I<0. Subsequently, the modulation of the parameter η₁(I) due to the decrease of I can be expected to induce a rate-induced resurgence of the AMOC. On the one hand, this is because the short sea ice timescale leads to very fast dynamics as the sea ice tips. On the other hand, even if the sea ice does not change quickly, when the amplitude of the change in η₁(I) becomes larger, there will be rate-induced tipping anyway due to the soft tipping point in the Stommel model described earlier. We thus choose the robust scenario where the coupling κ is such that the ocean component remains in the bi-stable regime with respect to η₁(I), and thus a rate-induced AMOC resurgence is the only pathway to tipping. As described in Sect. 2.2, this can be exemplified by a change in η₁(I) from η₁=2.65 (at the stadial sea ice fixed point for R=0) to η₁=3.0 for a collapsed sea ice cover I<0. Simulations with these parameters are qualitatively representative for a wider range of coupling strengths and rates of changing R.

Figure 10 Cascading stadial–interstadial transition in the coupled sea ice–ocean model where R is ramped from R=0 to $R=-\mathrm{0.3}$ within 340 years and kept constant afterwards. (a) Trajectory of the sea ice component as function of the control parameter R. (b, c) Trajectory of the ocean component as function of the changing parameter ${\mathit{\eta }}_{\mathrm{1}}\left(I\right(t\left)\right)={\mathit{\eta }}_{\mathrm{1}}^{\mathrm{0}}-\mathit{\kappa }\mathrm{\Theta }\left(I\right(t\left)\right)\cdot I\left(t\right)$. The tipping cascade consists of several steps separated by the points A, B, and C and marked by different colors in the trajectories (see main text). The gray surface in (c) is the moving basin boundary corresponding to the changing η₁(I(t)).

Figure 11 Probability of a cascading transition in the coupled sea ice–ocean model when changing the control parameter R linearly from R=0 to $R=-\mathrm{0.3}$ within different ramping times. (a) Fixed noise level σ_I=0.02 in the sea ice component and varying noise levels ${\mathit{\sigma }}_T={\mathit{\sigma }}_S=\mathit{\sigma }=\mathrm{0.005}$ (lightest gray), σ=0.01, σ=0.02, σ=0.04 and σ=0.1 (darkest gray) in the ocean component. (b) Fixed noise level ${\mathit{\sigma }}_T={\mathit{\sigma }}_S=\mathrm{0.02}$ in the ocean component and varying noise levels σ_I=0.005 (lightest gray), σ_I=0.01, σ_I=0.02, σ_I=0.04, σ_I=0.08 and σ_I=0.2 (darkest gray) in the sea ice component.

Figure 10 shows trajectories for a cascading stadial–interstadial transition in the deterministic limit when R is ramped down from R=0 to $R=-\mathrm{0.3}$ over 340 years. The transition can be divided into several stages. First, the sea ice slowly decreases as R is decreased and the ocean component tries to track the moving equilibrium (green segment of trajectories in Fig. 10). At point A, 325 years after the start of ramping, the sea ice passes the bifurcation point and rapidly tips to I<0 (purple segment in Fig. 10). This leads to a quick movement of η₁(I) towards η₁=3.0, which is reached at point B, 350 years after the start of ramping (Fig. 10b). As a result, the ocean state crosses the moving basin boundary (gray surface in Fig. 10c) from above and is thus determined to undergo rate-induced tipping to the on attractor (solid black curve). However, before tipping the ocean state is attracted by the stable manifold (i.e., the basin boundary) of the saddle (yellow segment). Finally, at point C (700 years after the start of ramping) the ocean component escapes the vicinity of the saddle and tips towards the on state (blue segment).

There is a critical timescale below which such a cascading transition with a rate-induced tipping is possible. This is a combination of the rate of change in the control parameter R and the speed of the tipping of the sea ice, which is held fixed here. As additive noise is included in the model, the boundary of tipping in terms of the ramping time of the control parameter is again blurred. Figure 11a shows the tipping probabilities for different noise levels as a function of the ramping time of R. The result is very similar to the ocean-only case, except that because of the fast tipping in the sea ice, the average ramping times leading to tipping are slightly higher. The picture looks different as we increase the noise level in the sea ice component, as seen in Fig. 11b. Here, the ramping times that yield significant tipping probabilities simply increase with the noise level without a large simultaneous decrease of the tipping probability for lower ramping durations. This is because noise-induced transitions to I<0 occur before the bifurcation of I is crossed. Since these transitions happen on the fast sea ice timescale, a rate-induced tipping of the ocean model becomes possible even when R is changed very slowly. As in the ocean-only case, the tipping cascade involves a saddle escape, which can lead to significant tipping delays as noise forcing of intermediate strength is included. Next, we will discuss this in more detail and relate it to potential pre-cursor signals leading up to such transitions.

3.4 Early warning of the tipping cascade

Due to their irreversible nature, it is important to foresee impending tipping points using generic early warning signals that do not require detailed knowledge of the system dynamics. These are typically obtained from time series by estimating a statistical indicator in a sliding window with appropriate detrending (see Sect. S4). For bifurcation tipping, a system often exhibits critical slowing down, which can be measured by increasing variance and autocorrelation. In Fig. 12 we show these indicators estimated in a sliding window for the cascading transition in Fig. 10. As expected there is an increase in variance and autocorrelation of I leading up to the bifurcation (Fig. 12c and d). Because of the speed of the parameter shift necessary to induce the cascade, the increases in the indicators do not fully exceed the variability prior to the parameter shift but could still provide early warning with a reasonable skill. Due to the coupling one might expect a signature of the sea ice critical slowing down in the ocean component. This is not seen here (Fig. 12e and f), since increasing fluctuations due to the sea ice are small compared to the variability in the ocean component for the chosen σ. If no noise is added to the ocean variables, critical slowing down can be detected in T or S. This might be an example of scenarios proposed in Rypdal (2016) and Boers (2018), where it is hypothesized that a bifurcation in the sea ice system is detectable as increased variance in the high frequencies of ice core data prior to DO events. Similarly, the fluctuations in I of increasing amplitude and temporal correlation may influence the ocean subsystem in a more consistent way as the bifurcation is approached. This should increase the cross-correlation, especially on longer timescales that can be measured with detrended cross-correlation analysis (DCCA). This has been proposed as early warning indicator for cascading transitions (Dekker et al., 2018). The method is similar to detrended fluctuation analysis, but instead of scaling in the variance, it measures scaling in the covariance of two signals with increasing timescales (for details, see Zebende, 2011 or Dekker et al., 2018). We can detect a slight increase on average in the DCCA exponent of I and T (Fig. 12g) for the transition in Fig. 10. However, the increase found in individual time series is not statistically significant, owing to the large variance of the DCCA estimator.

Figure 12 Ensemble simulations of the coupled sea ice ocean model, where R is ramped linearly from R=0 to $R=-\mathrm{0.3}$ within 350 years. (a) Time series of R (dashed line) and mean time series of I with a 90 % confidence band of the ensemble (gray shading). (b) Mean time series and 90 % confidence band of T and S. (c–f) Indicators of critical slowing down for I and T, estimated in a sliding window of 150 years, where the data in the window is detrended by a cubic function. The data is cut as the bifurcation in I is crossed until after the last realization tips plus the sliding window length. (g) Detrended cross-correlation analysis (DCCA) exponent estimated from I and T.

Figure 13 (a) Simulation in phase space of the Stommel model with ${\mathit{\sigma }}_T={\mathit{\sigma }}_S=\mathrm{0.2}$, where η₁ is ramped from η₁=2.65 to η₁=3.00 within 300 years. The two dotted lines correspond to the levels $q=T-S=\mathrm{0.06}$ and $q=T-S=\mathrm{0.1}$. The trajectory in between the first crossing of these two thresholds is shown in yellow. (b) Corresponding time series of the variable T.

3.5 Early warning of rate-induced tipping in the Stommel model

During the rate-induced transition of the ocean component there is an increase in the ensemble variance, as can be seen by the shadings in Fig. 12b. This increase, as well as a corresponding increase in ensemble autocorrelation, has been proposed before as an early warning signal for rate-induced tipping (Ritchie and Sieber, 2016). However, we show here that this results from the large spread in the amount of time spent by individual realizations at the saddle before tipping to the other attractor (see Fig. 8). In contrast, the fluctuations in individual realizations, as used for operational early warning, do not show an increase in variance and autocorrelation. This can be seen in Fig. 12e and f, where no increases in sliding window variance and autocorrelation accompany the increase in ensemble variance. For the estimation of variance and autocorrelation in a sliding window, a detrending of the time series is necessary, such that remaining trends in the residuals are not larger than the fluctuations themselves. For our detrending method using cubic functions, the severity of detrending, and thus the ability to remove sharp changes in the signal trend, depends only on the sliding window size (see Sect. S4 for more details). In order to remove the trend due to the parameter shift regarded here, a window size of no more than 200 years is required (Fig. S4).

Detrending inevitably removes some of the original fluctuations. To show that the lack of increased fluctuations in the detrended time series is not a consequence of too severe detrending, we extract segments of time series simulated from the Stommel-only model, where the system is in the vicinity of the saddle and there are no sharp trends. The fluctuations around the saddle are then compared to time series segments where the system fluctuates around the initial attractor. We define the vicinity of the saddle by the time periods in the simulations where q=0.06 is first crossed until q=0.1 is first crossed (Fig. 13). We regard the time series segments of an ensemble of realizations where the system stays in this vicinity for at least a certain duration. After detrending the segments by cubic functions, we calculate variance and autocorrelation. This yields empirical distributions of these quantities, describing the fluctuations in the system shortly before tipping. For each realization, we also choose a segment of the same duration taken just before the parameter shift starts, yielding distributions of variance and autocorrelation at the initial attractor. Figure 14 shows that variance and autocorrelation at the saddle are not increased but actually slightly decreased compared to the initial attractor. This is best seen for longer segments (Fig. 14c and d), since there the uncertainty in the estimators is smaller. One can also see that in this case the average variance and autocorrelation are larger compared to Fig. 14a and b because the detrending in longer windows removes less variability on longer timescales.

Figure 14 Distributions of variance and autocorrelation for ensembles of time series from the Stommel model (${\mathit{\sigma }}_T={\mathit{\sigma }}_S=\mathrm{0.2}$). These are estimated from time series segments around the initial fixed point at η₁=2.65 (black) and close to the saddle point (orange; see main text). The length of the segments for each realization corresponds to the time period that the system spent in the vicinity of the saddle. (a, b) Results for realizations where these time windows were at least 300 years long. (c, d) Results for time windows of at least 700 years. Also shown are the distributions around the on attractor (dashed) and the off attractor at η₁=3.0 (dotted).

Figure 15 (a–f) Distributions of the early warning indicator 𝒥 for ensembles of time series from the Stommel model (${\mathit{\sigma }}_T={\mathit{\sigma }}_S=\mathrm{0.2}$), estimated around the initial fixed point at η₁=2.65 (black) and close to the saddle point (orange). For each realization, 𝒥 is estimated after detrending in a time window that corresponds to the time period that the system spent in the vicinity of the saddle. In increasing order, the panels show results for realizations where these time windows were at least 100, 150, 200, 300, 400, and 600 years long, respectively. (g) Receiver operator characteristic curves for the same time series ensembles, showing the false and true positive rates as the threshold 𝒥_c is increased from low (top right) to high values (bottom left). The increasing darkness in the gray scale of the curves corresponds to the increasing time window lengths, as above. The diagonal dashed line indicates the performance of a pure chance classifier. The red cross indicates a perfect classifier.

It thus does not appear that critical slowing down indicators apply to rate-induced tipping. Instead, we exploit that the system is attracted towards the saddle where the dynamics are different to those at the initial attractor. If this difference can be detected before the system tips, a small perturbation in the right direction or a reversal of the parameter shift could push the system back in the desired basin of attraction. Saddles, which have at least one unstable direction in phase space, can be distinguished from attractors by a change from a negative to a positive real part of the largest eigenvalue of the Jacobian. Estimating the Jacobian from the time series in a sliding window could thus be a generic tool to detect the saddle escape involved in rate-induced tipping, and we describe a method to do this in the Appendix A. With this method the elements of the Jacobian during rate-induced tipping of the Stommel model can be inferred and allow for the distinction of the dynamics around the different fixed points (Fig. S5). However, there are quantitative biases in the estimates of individual elements, and as a result the estimates of the real part of the largest eigenvalue in the vicinity of the saddle are not consistently positive. These biases could be a result of the detrending, of a too high noise level or because the unstable dynamics are “suppressed” since we consider time series segments taken before the escape from the saddle.

As a more reliable indicator we propose the actual elements of the Jacobian, since they are inferred in a qualitatively robust way (Fig. S5). This lowers the estimator variance compared to the eigenvalues, which are composed of the estimates of all elements. The off-diagonal elements record changes in sign of the feedbacks in between the system variables. Such changes in feedback are common as a system moves towards a saddle. We combine the off-diagonal elements to a scalar early warning indicator 𝒥, as defined in Eq. (A6). Figure 15a–f shows that 𝒥 can distinguish the dynamics around the attractor (black) and the saddle (red) before tipping. The panels correspond to different minimum lengths of the time windows used to estimate 𝒥. The figure also shows probabilities p of observing a value of 𝒥 estimated around the attractor that is larger than a value of 𝒥 in the vicinity of the saddle. This measures the performance of 𝒥 as an early warning signal. For longer time windows, the distributions become better separated since the uncertainty of the estimator is reduced. Still, even for relatively short windows the indicator correctly identifies the departure from the attractor for most realizations.

Figure 16 Early warning indicators estimated in a 200-year sliding window from an ensemble of time series of the Stommel model, where η₁ is ramped from η₁=2.65 to η₁=3.0 within 300 years. (a) Time series of T and the parameter ramp. (b) Variance estimated from the detrended time series, as well as the ensemble variance (orange). (c) Lag 1 autocorrelation in the sliding window. (d) Early warning indicator 𝒥 (Eq. A6) estimated from the Jacobian in the sliding window. Mean time series are shown in black and the range in between 5th and 95th percentiles are shaded in gray.

An operational early warning signal can be constructed by estimating 𝒥 in a sliding window, and raising an alert as soon as a threshold 𝒥_c is exceeded. Choosing a location of 𝒥_c relative to the tails of the distributions in Fig. 15a–f is a trade-off in between maximizing the rate of true positives and minimizing the rate of false positives (alerts). The performance of the alert as a binary classifier can be summarized in receiver operating characteristic (ROC) curves. The curve of a perfect classifier collapses to the point (0,1). Figure 15g shows that for realizations that spend a longer time at the saddle, the indicator 𝒥 comes close to a perfect classifier, detecting the saddle approach with very low false positive and very high true positive rates. Figure 16 shows 𝒥 estimated from time series in a sliding window, along with critical slowing down indicators. 𝒥 begins to rise sharply roughly 200 years after the ramping started and decreases slightly as most realizations leave the saddle towards the on attractor. In contrast to the ensemble variance (orange), the variance and autocorrelation in the sliding window show no signal, apart from a small artifact around the parameter shift, which is a remnant of imperfect detrending in the 200-year windows.

4 Discussion

In this work we propose a conceptual model that describes a mechanism for abrupt climate change comprising a rate-induced resurgence of the AMOC as a response to increasing atmosphere–ocean heat exchange, which results from fast disappearance of sea ice. The latter occurs via a bifurcation tipping as a response to changing sea ice export into the North Atlantic, which could be driven by changes in wind stress forcing due to variations in ice sheet topography. In the context of DO events, the proposed model merely describes the sequence of events leading to a stadial–interstadial transition and not the dynamics of entire DO cycles that repeat in a self-sustained way. The model omits processes on longer timescales, as well as processes that would initiate after the resurgence of the AMOC. However, it can be easily extended to display self-sustained DO cycles by adding another slow variable that dynamically models the parameter shift. This could be a simple negative feedback reflecting, e.g., the influence of the AMOC on the ice sheets. Similarly, stronger noise forcing of the sea ice together with a weak feedback from the ocean to the sea ice can yield an excitable system with stochastically driven DO cycles. The proposed mechanism is thus a dynamical skeleton that is in principle compatible with stochastic, externally forced and self-sustained oscillatory dynamics driving DO cycles. Whether it indeed played a role in past abrupt climate change remains to be confirmed with more complex models, as well as with analyses of new highly resolved and synchronized climate proxy records.

The type of cascade introduced here could be a common feature in coupled systems that feature multistability and timescale separation. Here, a tipping in a fast subsystem can trigger a rate-induced transition of a slower subsystem even for weak coupling. Conversely, when there is no timescale separation but stronger coupling, the cascade can still occur in systems with non-smooth fold bifurcations. This is due to soft tipping points (Sect. 3.1), where the critical ramping duration to enable rate-induced tipping diverges as the parameter shift increases towards the bifurcation point. As a result, the cascading dynamics seen in the conceptual model may also be relevant for other regime shifts in the climate system, as well as for other natural systems. Consequently, we examined the mathematical details of the tipping cascade, which occurs in several stages. During the parameter shift the ocean subsystem tries to track the moving equilibrium. As the sea ice component tips abruptly, this fails and the system is instead attracted by the stable manifold of the saddle. The system then remains in the vicinity of the saddle as the dynamics slow down, before escaping to the on attractor. Adding noise leads to a broad distribution of the tipping time towards the on attractor. Early tipping, where stochastic perturbations push the system away from the stable manifold, is observed as well as significantly delayed tipping. In the latter case, noise pushes the system very close to the saddle, where it can get stuck for a very long time. A similar delay of rate-induced tipping for low noise levels has been reported for a one-dimensional gradient system (Ritchie and Sieber, 2016). It is seen from our model that due to the attraction by the stable manifold of the saddle, the tipping delay is a robust feature that exists for a fairly large range of rates (both sub- and super-critical) and noise levels. Thus, it opens up the possibility of issuing an early warning for rate-induced transitions.

The main difficulty for achieving an early warning of the cascade before the initial tipping of the sea ice is due to the relatively fast parameter shift involved. Thus, indicators proposed for cascading tipping points (Dekker et al., 2018) yield non-significant results, and more research is needed to find better indicators that might rely on similar principles. Instead, we focused on the rate-induced tipping of the ocean subsystem, since early warning signals for rate-induced tipping have not been developed. As in the case of fast passages through a bifurcation, for very fast parameter shifts one cannot hope for an early warning of rate-induced tipping. Here the system is not attracted by the saddle but evolves quickly towards the alternative attractor. However, for intermediate rates we can exploit the fact that tipping occurs via saddle escape. As the system state departs from the moving attractor towards the saddle, the linear stability changes. This can be captured by the Jacobian matrix, which we estimate from the time series. We then propose using the off-diagonal elements of the Jacobian as an early warning signal. These elements record changes in the sign of coupling in between the system variables, indicating a change in stability. The proposed indicator detects an approach of the saddle with significant skill, in particular for realizations where the system stays in the vicinity for a longer time, so that the Jacobian can be estimated with good precision. Note that the actual tipping occurs by escaping the vicinity of the saddle, which is largely noise-induced. Thus, early warning in the sense of predicting the precise time of the saddle escape is hard to achieve. Early-warning signals for saddle escapes have been proposed (Kuehn et al., 2015), but they require being very close to the saddle and very low noise. While the specific early warning signal proposed here may not apply to all cases of rate-induced tipping, the general procedure of detecting a qualitative change in the feedback structure of the system via the Jacobian or its eigenvalues should be widely applicable. For higher-dimensional systems early warning might even become easier, since there are often dominant eigenvalues and large differences in the effective dimensionality of the dynamics on the attractor vs. the transient dynamics during tipping. Other techniques for detecting transient dynamics may also be useful here (Gottwald and Gugole, 2020). The phenomenology of cascading transitions involving rate-induced tipping that has been exemplified here is to be tested with models of different complexity in upcoming studies. Furthermore, the applicability of the early warning method to real-world data needs to be tested. In the typical case where only one (or a few) scalar time series are available, this will involve a time series embedding and subsequent estimation of the Jacobian from the reconstructed multivariate time series.

5 Conclusions

Building on previous studies of proxy records and state-of-the-art climate models, we propose that past abrupt climate change could have arisen as a cascade of tipping points. We translate this into a conceptual sea ice–ocean model, where a parameter shift leads to the following cascade. First, as a result of the gradually changing climatic conditions, the North Atlantic sea ice cover collapses abruptly. Subsequently, the AMOC resurges abruptly from a weak to a vigorous state in a rate-induced tipping, as a response to the fast rate of sea ice decline enhancing the atmosphere–ocean heat exchange. Our analysis suggests that cascades of tipping points in weakly coupled climate components with timescale separation become more likely under certain circumstances. This is the case when there are rate-dependent tipping points or soft tipping points associated with non-smooth fold bifurcations. This motivates the development of specialized early warning signals for such rate-dependent cascading tipping points. We present a first step in this direction by showing that due to a delay in the tipping of the ocean circulation a statistical estimation of the Jacobian can detect the impending abrupt transition. This may be applicable as generic early warning signal of rate-induced transitions.

Appendix A: An early warning indicator for rate-induced tipping

We detect rate-induced tipping by identifying a departure from the initial attractor towards the vicinity of the saddle. This is accompanied by a change in the linear stability of the system and thus the Jacobian. The latter is estimated from the multi-variate time series in a sliding window as follows. Consider the underlying dynamical system $\dot{\mathit{x}}\left(t\right)=\mathit{f}\left(\mathit{x}\right(t\left)\right)$ with x∈ℝ^d and the observed discrete time series $\mathit{\left\{}\mathit{x}\right(\mathrm{1}),\mathit{x}(\mathrm{2}),\mathrm{\dots },\mathit{x}(N\left)\mathit{\right\}}$, where N is the window size. The linearization of the dynamical system around the equilibrium point y is

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\displaystyle }\dot{\stackrel{\mathrm{̃}}{\mathit{x}}}\left(t\right)=\sum _{i=\mathrm{1}}^d\stackrel{\mathrm{̃}}{x_i}\left(t\right){\displaystyle \frac{\partial \mathit{f}\left(\mathit{x}\right)}{\partial x_i}}{|}_{\mathit{x}=\mathit{y}},\end{array}$$

with $\stackrel{\mathrm{̃}}{\mathit{x}}\left(t\right)=\mathit{x}\left(t\right)-\mathit{y}$. Discretized, this can be approximated as

$$\begin{array}{ll}{\displaystyle }\stackrel{\mathrm{̃}}{\mathit{x}}(t+\mathrm{1})-\stackrel{\mathrm{̃}}{\mathit{x}}\left(t\right)& {\displaystyle }=\mathit{x}(t+\mathrm{1})-\mathit{x}\left(t\right)\equiv \mathrm{\Delta }{\mathit{x}}_t\\ \text{(A2)}& {\displaystyle }& {\displaystyle }=\mathit{\delta }t\left(\sum _{i=\mathrm{1}}^d\stackrel{\mathrm{̃}}{x_i}\left(t\right){\displaystyle \frac{\partial \mathit{f}\left(\mathit{x}\right)}{\partial x_i}}{|}_{\mathit{x}=\mathit{y}}\right).\end{array}$$

In this expression, the factors $\frac{\partial f_j\left(\mathit{y}\right)}{\partial x_i}$ are the elements J_ji of the Jacobian matrix. They can be estimated with multiple linear regression by sampling different Δx_t as a dependent variable and $\stackrel{\mathrm{̃}}{x_i}\left(t\right)$ for i=1 … d as independent variables for a given y from within the time series. To this end, we choose $\mathit{\left\{}\mathit{x}\right(t_{\mathrm{1}}),\mathit{x}(t_{\mathrm{2}}),\mathrm{\dots },\mathit{x}(t_M\left)\mathit{\right\}}$ from within the windowed time series, which are the M closest points to y in phase space in terms of the distance $D_{y,k}={\sum }_{i=\mathrm{1}}^d\left[x_i\right(t_k)-y_i{]}^{\mathrm{2}}$. For each x(t_k), we evaluate $\mathrm{\Delta }{\mathit{x}}_{t_k}$ using the subsequent point in the time series. From the M samples of $\mathrm{\Delta }{\mathit{x}}_{t_k}$ and $\stackrel{\mathrm{̃}}{x_i}\left(t_k\right)$ for i=1…d, we obtain the factors $\frac{\partial f_j\left(\mathit{y}\right)}{\partial x_i}$ by multiple linear regression. We then repeat the procedure for every data point in the window as y, and average the results to obtain average Jacobian elements J_ji within the sliding window. In this work we chose $M=N/\mathrm{2}$. To illustrate how the Jacobian changes in the Stommel model as the system departs the off attractor, we write Eq. (4) in the deterministic case as

$$\begin{array}{}\text{(A3)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=f(T,S),\text{(A4)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=g(T,S).\end{array}$$

The corresponding Jacobian of the linearized system is as follows.

$$\begin{array}{}\text{(A5)}& \begin{array}{l}\mathbf{J}=\left(\begin{array}{cc}\frac{\partial f(T,S)}{\partial T}& \frac{\partial f(T,S)}{\partial S}\\ \frac{\partial g(T,S)}{\partial T}& \frac{\partial g(T,S)}{\partial T}\end{array}\right)\\ =\left(\begin{array}{cc}\text{sgn}(T-S)\cdot (S-\mathrm{2}T)-\mathrm{1}& \text{sgn}(T-S)\cdot T\\ \text{sgn}(S-T)\cdot S& \text{sgn}(T-S)\cdot (\mathrm{2}S-T)-{\mathit{\eta }}_{\mathrm{3}}.\end{array}\right)\end{array}\end{array}$$

Around the attractors, the real parts of both eigenvalues are negative. As the saddle is approached by crossing q>0, the real part of the first eigenvalue becomes positive. Furthermore, the off-diagonal elements of the Jacobian change sign. We propose this sign change as early warning signal, since it is more robust than the eigenvalues when estimated from noisy data. We define the early warning signal as

$$\begin{array}{}\text{(A6)}& {\displaystyle }{\displaystyle }\mathcal{J}\equiv {\displaystyle \frac{\partial f}{\partial S}}-{\displaystyle \frac{\partial g}{\partial T}}.\end{array}$$

Note that for dynamical systems defined by a gradient of a potential this indicator is not applicable, since it would be 0 in the whole phase space due to the symmetric Jacobian. Using just one of the diagonal elements as indicator instead still gives good early warning possibilities with roughly half the statistical power due to the smaller amount of information retained. For time series from unknown dynamical systems, changes in the individual elements could be monitored simultaneously, potentially after embedding in the case of univariate time series.