Cascading transitions in the climate system

We introduce a framework of cascading tipping, i.e. a sequence of abrupt transitions occurring because a transition in one subsystem changes the background conditions for another subsystem. A mathematical framework of elementary deterministic cascading tipping points in autonomous dynamical systems is presented containing the double-fold, fold–Hopf, Hopf–fold and double-Hopf as the most generic cases. Statistical indicators which can be used as early warning indicators of cascading tipping events in stochastic, non-stationary systems are suggested. The concept of cascading tipping is illustrated through a conceptual model of the coupled North Atlantic Ocean – El Niño–Southern Oscillation (ENSO) system, demonstrating the possibility of such cascading events in the climate system.

1 Introduction

Earth's climate system consists of several subsystems, e.g. the ocean, the atmosphere, ice and land, which are coupled through fluxes of momentum, mass and heat. Each of these subsystems is characterized by specific processes, on very different timescales, determining the evolution of its observables. For example, processes in the atmosphere occur on much smaller timescales than in the ocean; hence, in weather prediction the upper ocean sets the background state for the evolution of the atmosphere. Similarly, in equatorial ocean–atmosphere dynamics associated with the El Niño–Southern Oscillation (ENSO) phenomenon, the global meridional overturning circulation can be considered a background state, as it evolves on a much larger timescale.

This notion that one subsystem provides a background state for the evolution of another subsystem is important when critical transitions are considered. In the climate system, a number of tipping elements have been identified (Lenton et al., 2008 give an overview of these), where changes in observables can occur relatively rapidly compared to the changes in their forcing near so-called “tipping points”. Examples of tipping elements are the Atlantic meridional overturning circulation (AMOC) (Stommel, 1961), the Arctic sea ice (Bathiany et al., 2016), monsoon patterns, mid-latitude atmospheric flow (Barriopedro et al., 2006), vegetation cover (Hirota et al., 2011) and more local systems like coral reefs and permafrost. When one subsystem undergoes a transition, which changes the background state of another subsystem, a transition may also be induced in that second subsystem. Such dynamical interactions which lead to coupled transitions are examples of “tipping cascades” or “domino effects” (Kriegler et al., 2009; Lenton and Williams, 2013).

Many tipping points have been analysed in separate subsystems, both for phenomena of the present-day climate (Bathiany et al., 2016; Lenton, 2011), and for past climates (such as the abrupt cooling of the Younger Dryas, Livina and Lenton, 2007, and the desertification of the Sahel region, Kutzbach et al., 1996). However, less attention has been paid to the interaction between transitions in different subsystems. For example, when the AMOC collapses, precipitation patterns may change such that the equilibrium structure of the vegetation cover in the Amazon rainforest is shifted (Aleina et al., 2013). This may result in another transition, concerned with forest growth or dieback. Another example is the influence of the AMOC on the trade winds (through meridional sea surface temperature gradients), which in turn influence the amplitude of ENSO. In models, a collapse of the AMOC has been found to intensify ENSO (Dong and Sutton, 2007; Lenton and Williams, 2013; Timmermann et al., 2007), although there are also other effects that would weaken ENSO (Timmermann et al., 2005).

An example in past climates is the coupling between the ocean's overturning circulation and land ice. The rapid glaciation of the Antarctic continent around the Eocene–Oligocene boundary (34 Ma) is often explained in terms of a CO₂ threshold being reached that allowed a major ice sheet to grow (DeConto and Pollard, 2003; Gasson et al., 2014). However, a two-step signal is found in the oxygen isotopic ratio, δ¹⁸O, which is attributed to a deep-sea temperature drop followed by the (slower) growth of the Antarctic Ice Sheet (AIS). One suggestion to explain the two-step transition is that the deep-sea temperature drop was related to a change in the pattern of the global MOC (Tigchelaar et al., 2011). The ice sheet formation is then argued to have been driven by decreasing atmospheric CO₂ (Pearson et al., 2009). The switch in MOC (first tipping) led to the changes in the atmospheric CO₂ (e.g. Elsworth et al., 2017) which caused the growth of the AIS (second tipping). This leads to the question of whether a cascading tipping event occurred.

In the last few years, much work has been carried out to formulate statistical indicators and early warning signals of tipping points. A system close to a critical transition shows features of a “critical slowing down” (Dakos et al., 2008; Kuehn, 2011; Scheffer et al., 2009). In the vicinity of the tipping point, the system slowly loses its ability to recover from small perturbations. This results in increased variance, auto-correlation and potentially also increased skewness and flickering (Scheffer et al., 2009). Various methods providing a specific scalar and a threshold when approaching the transition have been suggested, such as degenerate fingerprinting (Held and Kleinen, 2004; Thompson and Sieber, 2011) and detrended fluctuation analysis (DFA) (Livina and Lenton, 2007; Peng et al., 1994).

When considering cascading tipping points, the auto-correlation of two time series and their interaction need to be analysed simultaneously. Podnobik and Stanley (2008) proposed an altered form of DFA to assess the cross-correlation between two non-stationary time series and called this method detrended cross-correlation analysis (DCCA). In the computation of the fluctuation function, they used cross-covariance instead of auto-covariance and fit this to a power law. This concept is further extended by defining a coefficient ρ_DCCA that accounts for the auto-covariance of the individual time series (Yuan et al., 2015; Zhou, 2008). However, no statistical analysis or indicators have yet been formulated for cascading tipping events.

In this paper, we provide a quantitative approach to cascading tipping events. We start with a mathematical framework to formulate elementary cascading tipping points (Sect. 2). Next, we discuss statistical metrics to analyse, diagnose and potentially predict cascading transitions, and apply them to ensemble simulations of the elementary cascading tipping points (Sect. 3). Finally, we apply the new concepts to an example within the climate system: the potential cascading tipping mechanism between the AMOC and ENSO (Sect. 4). We summarize and discuss our findings in Sect. 5.

2 Mathematical framework for cascading tipping

In the climate system, tipping points are usually related to rapid transitions, where an observable in the climate system may change abruptly in a relatively short time compared with changes in the forcing of the observable. Such rapid changes usually involve transitions from one equilibrium state to another, which can often be explained using classical bifurcation theory for autonomous dynamical systems. To a certain extent, these concepts can also be applied to non-autonomous systems (so-called slow–fast systems) when the time variation of parameters can be viewed as a slow external forcing (Kuehn, 2011). They form the basics for understanding phenomena such as noise-induced tipping (Thompson and Sieber, 2011) and rate-dependent tipping (Ashwin et al., 2012).

In this section, we present a mathematical framework for simple cascading transitions, which acts as a first step towards analysing the more complex real-world transitions. We focus on bifurcation-induced tipping points, and consider two types of bifurcations that are thought to be relevant to mechanisms of abrupt changes in the climate system; the back-to-back saddle-node bifurcation is often used to explain transitions between two coexisting equilibria (multi-stable systems), while the Hopf bifurcation can explain the appearance of oscillatory behaviour (Thompson and Stewart, 2002). Therefore, abrupt change in the system appears as a consequence of a parameter crossing a specific critical value at the bifurcation point.

A back-to-back saddle-node bifurcation (two saddle-nodes connected by a common unstable branch) generically occurs in physical systems (that have bounded states) when one parameter varies and the simplest dynamical system has a bifurcation that is described by

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}=a_{\mathrm{1}}{x}^{\mathrm{3}}+a_{\mathrm{2}}x+\mathit{\varphi },\end{array}$$

where the a_i ($i\in \mathrm{1},\mathrm{2}$) are constants, ϕ is a parameter, x is the state variable and t is time. There are multiple equilibria in the system if and only if a₁<0, a₂>0 and within the parameter interval $\left|\mathit{\varphi }\right|<\left(\right(-\mathrm{4}{a}_{\mathrm{1}}^{\mathrm{3}}{a}_{\mathrm{2}}^{\mathrm{3}})/(\mathrm{27}{a}_{\mathrm{1}}^{\mathrm{4}}){)}^{\mathrm{1}/\mathrm{2}}$. In this case, the back-to-back saddle-node bifurcations occur at ${\mathit{\varphi }}_c=\pm \left(\right(-\mathrm{4}{a}_{\mathrm{1}}^{\mathrm{3}}{a}_{\mathrm{2}}^{\mathrm{3}})/(\mathrm{27}{a}_{\mathrm{1}}^{\mathrm{4}}){)}^{\mathrm{1}/\mathrm{2}}$. In the following, we often use the terms “fold bifurcation” and “saddle-node bifurcation”, although in practise we only use “back-to-back” saddle-node bifurcations.

Table 1 Parameter values and coupling for the four types of cascading tipping as shown in Figs. 1 and 2.

A Hopf bifurcation also generically occurs in physical systems, and the simplest dynamical system in which it occurs when one parameter is varied is described by

$$\begin{array}{}\text{(2)}& {\displaystyle }\left\{\begin{array}{l}\frac{\mathrm{d}x}{\mathrm{d}t}=a_{\mathrm{1}}y+a_{\mathrm{2}}(\mathit{\varphi }-(x^{\mathrm{2}}+y^{\mathrm{2}}\left)\right)x\\ \frac{\mathrm{d}y}{\mathrm{d}t}=b_{\mathrm{1}}x+b_{\mathrm{2}}(\mathit{\varphi }-(x^{\mathrm{2}}+y^{\mathrm{2}}\left)\right)y,\end{array}\right.\end{array}$$

where, again, $a_i,b_i,i=\mathrm{1},\mathrm{2}$ are constants, ϕ is the parameter, (x, y) is the state vector and t is time. The state vector satisfying (Eq. 2) reaches a stable periodic orbit if and only if a₁b₁<0 and ϕ>0; the transition from steady to periodic occurs at ϕ=0.

There are two other bifurcations when one parameter is varied (the transcritical and pitchfork bifurcations); however, these bifurcations are non-generic because special conditions must hold (e.g. symmetry) and so they are not considered here. Using the saddle-node and Hopf bifurcations, cascading tipping can be viewed as a combination of two coupled subsystems, where each subsystem undergoes one of these two types of bifurcations. Combining these bifurcations leads to four types of cascading tipping, which are discussed in the following.

Figure 1 Stable (solid), unstable (black dotted) and oscillatory (amplitude red dotted) regimes of various cascading tipping types (as defined on top of the figure). Red, orange and black dots indicate fold, Hopf and torus bifurcations, respectively. Leading system versus forcing ϕ (a–d). Following system versus coupling γ (e–h). Following system versus forcing ϕ (via leading) (i–l). Critical values of (coming from lower branch) leading tipping (a–d), following tipping (e–h) and the combined cascading tipping event (i–l) are marked by the grey vertical dashed line.

Figure 2 Example simulations for each cascading event type: the double-fold cascade (a), the fold–Hopf cascade (b), the Hopf–fold cascade (c) and the double Hopf cascade (d). Black and grey lines indicate the leading systems, red and yellow lines indicate the following systems. Dotted lines indicate time before the critical threshold in the forcing ϕ(t) (black/grey) or coupling κ(x) (red/yellow) is reached, solid lines indicate the time after this. Parameter values for the modelled systems are given in Table 1.

Coupling two systems introduces a direction to the cascade and we take account of this by defining a “leading” system, which during its transition changes a parameter (that is, the coupling term) in the “following” system. The changing parameter in the following system can then bring the following system closer to a bifurcation point, potentially even resulting in a second transition. In this section, we only look at deterministic cases, which do not allow for noise to play a role in the tipping. Therefore, transitions in the leading system result in a changed coupling term that can lead to transitions in the following system. In bifurcation diagrams versus forcing, the bifurcation points (for deterministic systems) can overlap. However, the transitions are distinguishable in transients, because the following system always tips after the completion of the first transition. This is why we show the bifurcation diagrams of both systems versus forcing (Fig. 1) and an example of a transient (Fig. 2) for each type. They will be discussed in more detail below.

2.1 Double-fold cascade

The most intuitive system that has the potential to undergo a cascading tipping event is a system where both the leading and the following system have saddle-node bifurcations. Analogous to Eq. (1), a dynamical system containing a double fold cascade is then

$$\begin{array}{}\text{(3)}& {\displaystyle }\left\{\begin{array}{l}\frac{\mathrm{d}x}{\mathrm{d}t}=a_{\mathrm{1}}{x}^{\mathrm{3}}+a_{\mathrm{2}}x+\mathit{\varphi }\\ \frac{\mathrm{d}y}{\mathrm{d}t}=b_{\mathrm{1}}{y}^{\mathrm{3}}+b_{\mathrm{2}}y+\mathit{\gamma }\left(x\right),\end{array}\right.\end{array}$$

where x is the state vector of the leading system, y is the state vector of the following system, a_i,b_i ($i\in \mathrm{1},\mathrm{2}$) are constants, and ϕ is a parameter in the leading system. The key here is that the function γ, which serves as a parameter in the following system, depends on the leading system. The most simple coupling between the two systems is represented by $\mathit{\gamma }\left(x\right)={\mathit{\gamma }}_{\mathrm{1}}+{\mathit{\gamma }}_{\mathrm{2}}x$. Observe that a change in the forcing parameter ϕ can induce a transition in x, which may affect the coupling function γ such that y also undergoes a transition. We would like to emphasize that the forcing ϕ does not directly affect y, and that y is only affected through a change in x (which is effectively only significant when x tips).

Implementing this system with the values shown in Table 1 gives insight in the system's equilibrium structure (Fig. 1a) and transient behaviour (example in Fig. 2a). When ϕ is changed moving through the bistable regime of the leading system the coupling moves the following system through its own bistable regime (see Table 1). Figure 1a and b show the equilibria of the leading (versus ϕ) and following systems (versus γ), respectively, showing the bistable regime in the centre of the figure, embedded in the saddle-node structure. Varying ϕ alters the state of the leading system, which affects the state of the following system through the coupling γ. This results in the existence of four stable states in the following system of the bistable regime of the leading system: two per state of the leading system, as shown in Fig. 1i. The leading system's state acts as a background condition modulating the position of the following system's equilibria; therefore, in case of transition, the leading system's state may drastically reposition the equilibria of the following system. This is intuitively visible in Fig. 2a, where a time series example of the dynamical system in Eq. (3) shows a cascading tipping event (parameters shown in Table 1). When the leading system (black) is forced (by changing ϕ) to move from a bistable to a monostable regime, it transits towards a new equilibrium. During this transition, the following system (red) is affected and leaves the regime in which it had four possible equilibria; it also transits to a different state.

2.2 Fold–Hopf cascade

The second type of cascading tipping event involves a saddle-node bifurcation in the leading system and a subsequent Hopf bifurcation in the following system. Using analogous notation as in Eq. (3), the simplest system that captures this so-called fold–Hopf cascade is

$$\begin{array}{}\text{(4)}& {\displaystyle }\left\{\begin{array}{l}\frac{\mathrm{d}x}{\mathrm{d}t}=a_{\mathrm{1}}{x}^{\mathrm{3}}+a_{\mathrm{2}}x+\mathit{\varphi }\\ \frac{\mathrm{d}y}{\mathrm{d}t}=b_{\mathrm{1}}z+b_{\mathrm{2}}\left(\mathit{\gamma }\right(x)-(y^{\mathrm{2}}+z^{\mathrm{2}}\left)\right)y\\ \frac{\mathrm{d}z}{\mathrm{d}t}=c_{\mathrm{1}}y+c_{\mathrm{2}}\left(\mathit{\gamma }\right(x)-(y^{\mathrm{2}}+z^{\mathrm{2}}\left)\right)z,\end{array}\right.\end{array}$$

where x is again the state vector of the leading system, and (y, z) is the state vector of the following system. By slowly varying the parameter ϕ (e.g. linearly as ϕ(t)) the leading system moves through its bistable regime (see Table 1 for parameter values) and via the coupling $\mathit{\gamma }\left(x\right)={\mathit{\gamma }}_{\mathrm{1}}+{\mathit{\gamma }}_{\mathrm{2}}x$ forces the following system across the Hopf bifurcation point.

The bifurcation structure of the leading system of Eq. (4), using the parameters stated in Table 1, is displayed in Fig. 1b. As in Fig. 1a, this system's bifurcation diagram shows a saddle-node structure. The following system in Fig. 1f, now undergoes a Hopf bifurcation and subsequent oscillatory behaviour. In Fig. 1j, it can be seen that increasing ϕ on the lower branch of the leading system regime will still result in a stationary equilibrium for the following system. When increasing ϕ across ϕ_c, the leading system moves towards another state (seen in Fig. 1b) and the coupling γ increases (across γ_c in Fig. 1f). An oscillation then occurs in the following system. This mechanism makes it possible for steady and oscillatory states to coexist on the right side of the Hopf bifurcation in Fig. 1j. An example of a time series showing a fold–Hopf cascading event is displayed in Fig. 2b. A transition in the leading system (black) brings the following system (red/yellow) into an unstable equilibrium that eventually leads to an oscillatory state.

2.3 Hopf–fold cascade

A third type of cascading event involves a Hopf bifurcation in the leading system and a subsequent saddle-node bifurcation in the following system. Using similar notation to the previous subsection, the simplest system with a Hopf–fold cascade (see Table 1 for parameter values) is given by

$$\begin{array}{}\text{(5)}& {\displaystyle }\left\{\begin{array}{l}\frac{\mathrm{d}x}{\mathrm{d}t}=a_{\mathrm{1}}y+a_{\mathrm{2}}(\mathit{\varphi }-(x^{\mathrm{2}}+y^{\mathrm{2}}\left)\right)x\\ \frac{\mathrm{d}y}{\mathrm{d}t}=b_{\mathrm{1}}x+b_{\mathrm{2}}(\mathit{\varphi }-(x^{\mathrm{2}}+y^{\mathrm{2}}\left)\right)y\\ \frac{\mathrm{d}z}{\mathrm{d}t}=c_{\mathrm{1}}{z}^{\mathrm{3}}+c_{\mathrm{2}}z+\mathit{\gamma }\left(x\right),\end{array}\right.\end{array}$$

where (x, y) is the state vector of the leading system, and z is the state vector of the following system. Again, we can slowly increase ϕ so that the leading system (x, y) crosses a Hopf bifurcation; via the coupling $\mathit{\gamma }\left(x\right)={\mathit{\gamma }}_{\mathrm{1}}+{\mathit{\gamma }}_{\mathrm{2}}x$ the following system is then moved through its bistable regime and a fold is reached in z.

Figure 1c contains the typical bifurcation structure of the leading system in Eq. (5), containing a Hopf bifurcation separating a stationary from an oscillatory regime. The following system's equilibrium structure for varying γ is given by Fig. 1g. In this particular configuration, for any negative ϕ there are multiple stable equilibria, as seen in Fig. 1k. This makes sense, as ϕ only affects the following system via its impact on the leading system, and for negative ϕ the leading system remains constant. At ϕ=0, the Hopf bifurcation in the leading system is reached and (x, y) starts oscillating. The following system oscillates a little along with the leading system due to the oscillatory changing value of γ.

When ϕ increases more, the amplitude of the leading system's oscillation grows, which eventually causes γ to cross the threshold; the following system then leaves its bistable regime (be it temporarily as γ will be reduced again due to the oscillation). This forces the following system into its upper branch, as can be seen in Fig. 1g. The upper branch's stable oscillation ends at ϕ≈0.8 (in Fig. 1k), because the amplitude becomes large enough for the system to swap between multiple equilibria. An example of such a cascading transition event can be seen in Fig. 2c, where an oscillation starts in the leading system (black/grey). A particular phase of this oscillation then causes the following system (red) to transit into the second equilibrium.

2.4 Double-Hopf cascade

A fourth type of cascading tipping event discussed here involves a Hopf bifurcation in the leading system and a subsequent Hopf bifurcation in the following system. Using analogous notation to the previous subsections, this double Hopf cascade is captured by the following dynamical system:

$$\begin{array}{}\text{(6)}& {\displaystyle }\left\{\begin{array}{l}\frac{\mathrm{d}x}{\mathrm{d}t}=a_{\mathrm{1}}y+a_{\mathrm{2}}(\mathit{\varphi }-(x^{\mathrm{2}}+y^{\mathrm{2}}\left)\right)x\\ \frac{\mathrm{d}y}{\mathrm{d}t}=b_{\mathrm{1}}x+b_{\mathrm{2}}(\mathit{\varphi }-(x^{\mathrm{2}}+y^{\mathrm{2}}\left)\right)y\\ \frac{\mathrm{d}u}{\mathrm{d}t}=c_{\mathrm{1}}v+c_{\mathrm{2}}\left(\mathit{\gamma }\right(x)-(u^{\mathrm{2}}+v^{\mathrm{2}}\left)\right)u\\ \frac{\mathrm{d}v}{\mathrm{d}t}=d_{\mathrm{1}}u+d_{\mathrm{2}}\left(\mathit{\gamma }\right(x)-(u^{\mathrm{2}}+v^{\mathrm{2}}\left)\right)v,\end{array}\right.\end{array}$$

where (x, y) is the state vector of the leading system, and (u, v) is the state vector of the following system. If ϕ forces (x, y) such that it crosses the Hopf bifurcation point, the coupling $\mathit{\gamma }\left(x\right)={\mathit{\gamma }}_{\mathrm{1}}+{\mathit{\gamma }}_{\mathrm{2}}x$ causes a crossing of the second Hopf bifurcation in (u, v).

Figure 1d and h show the bifurcation diagrams of the leading and following systems, with supercritical Hopf bifurcations. The following system (in Fig. 1l) is stationary for many values of ϕ, up to the point that the leading system starts oscillating, which for high enough values of ϕ is large enough to make the following system cross the Hopf bifurcation. However, γ oscillates with the leading system (for ϕ>0). This means that oscillatory behaviour in the following system can only be expected in a particular part of the leading system's oscillation period. This interaction between the two oscillations result in torus bifurcations for particular values of ϕ. An example of a time series showing a double-Hopf cascading transition is shown in Fig. 2d. After a (slow) oscillation in the leading system (black/grey) has started, a (fast) oscillation in the following system (red/yellow) arises in particular phases of the slow oscillation.

3 Statistics of cascading tipping points

In the previous section we formulated elementary deterministic dynamical systems that can exhibit cascading tipping. In order to detect tipping events from, for example, observed time series in real systems, we need to detect whether a system is close to a critical transition. In general, systems close to critical transition recover more slowly from perturbations, which in turn increases memory in the time series. This leads to the phenomenon of “critical slowing down” prior to bifurcation points. In this section, we simulate cascading tipping events to (a) study the statistical character of such events, (b) diagnose (post-tipping) whether tipping events are causally related and (c) take a first step towards statistical indicators for the prediction of cascading tipping events.

3.1 Methods for single tipping points

Several methods have been suggested for the analysis of time series to detect the approach of a single tipping point. For saddle-node bifurcations, the key features of a time series such as this is a critical slowing down. This can be investigated using standard quantities such as increasing auto-correlation (e.g. the lag-1 auto-correlation), increasing variance and increasing skewness (Held and Kleinen, 2004; Kuehn, 2011; Scheffer et al., 2009). Although critical slowing down near critical transitions is a general feature of (even chaotic) dynamical systems (Tantet et al., 2018), the standard metrics like auto-correlation at lag 1 and variance do not always provide an early warning signal (e.g. in Greenland ice core data in Livina and Lenton, 2007). Hence, more complicated indicators have been introduced, such as (i)  degenerate fingerprinting (DF) and (ii) the detrended fluctuation analysis (DFA) (Held and Kleinen, 2004; Livina and Lenton, 2007; Peng et al., 1994; Thompson and Sieber, 2011). DFA is argued to be a solution to the problem that the simple lag-1 auto-correlation does not capture the approach to a transition in highly non-stationary data in long time series (Livina and Lenton, 2007; Peng et al., 1994). The latter is characterized by high auto-correlation due to the gradual increase or decrease of the system itself, distorting the signal of the critical slowing down, a problem in standard metrics and DF.

As critical slowing down implies an increasing auto-regressive behaviour in the time series prior to a transition, the memory component is increased. In general, the first step in DF is the projection of the data fields onto the leading EOF, which gives a time series x_n (Held and Kleinen, 2004). After time-equidistant interpolation and detrending of the data, in the DF method one fits the following general auto-regressive process to the time series x_n:

$$\begin{array}{}\text{(7)}& x_{n+\mathrm{1}}=c\cdot x_n+\mathit{\sigma }{\mathit{\eta }}_n,\end{array}$$

where η_n is Gaussian white noise and $c=\text{exp}(-\mathit{\lambda }\mathrm{\Delta }t)$ is the AR(1) coefficient. Here λ can be seen as the decay rate of perturbations in previous time steps. As the approaching of a bifurcation point involves an increase in memory, the value of c is presumed to increase towards one when approaching a saddle-node bifurcation point.

In DFA, one first chooses an integer window size s and divides the (cumulative-summed) time series X(n) into $N_s=N/s$ segments that do not overlap, where N is the length of the time series. In every window, the best polynomial fit of a chosen order is calculated. A quadratic polynomial is used here. The squared deviation from this quadratic polynomial for every window is summed, resulting in a measure of the auto-covariance fluctuating around the fit:

$$\begin{array}{}\text{(8)}& {\displaystyle }{F}^{\mathrm{2}}(\mathit{\nu },s)={\displaystyle \frac{\mathrm{1}}{s}}\sum _{i=\mathrm{1}}^s\left[X\right((\mathit{\nu }-\mathrm{1})s+i)-x_{\mathit{\nu }}(i){]}^{\mathrm{2}},\end{array}$$

where X is the detrended time series and x_ν is the best polynomial fit in segment ν. Next, an average is taken over all segments to obtain the fluctuation function F(s):

$$\begin{array}{}\text{(9)}& {\displaystyle }F\left(s\right)=\sqrt{{\displaystyle \frac{\mathrm{1}}{N/s}}\sum _{\mathit{\nu }=\mathrm{1}}^{N/s}{F}^{\mathrm{2}}(\mathit{\nu },s)}\end{array}$$

which depends solely on s. The long-range auto-correlations can now be recognized by fitting the fluctuation function to a power law and looking at the resulting DFA-exponent α, according to

$$\begin{array}{}\text{(10)}& {\displaystyle }F\left(s\right)\propto s^{\mathit{\alpha }}.\end{array}$$

For α≤0.5, there is no long-term correlation and the fluctuations are indistinguishable from white noise. However, when α>0.5, there are long-term correlations present and for α≥1.5 the system has reached a bifurcation point. In the simulations analysed here the DFA scaling exponent is fitted explicitly for every (moving) window.

3.2 Methods for cascading tipping points

Cascading tipping involves two systems with their own bifurcation structure and their proximity towards bifurcation points. Although the leading system may be close to tipping, the following system might still be far from its bifurcation point and needs the critical transition of the leading system to even come close to this point; this makes the behaviour of the following system more prone to noise, less dependent on the leading system and less auto-correlated prior to the first tipping. This is why the general measures for single tipping events cannot be used, nor can regular cross (Pearson) correlation. The reason for this is that the following and leading systems do not have a one-to-one relationship (that is, weakly Pearson correlated), but are rather coupled through specific parameters, which is only visible in long-range correlations.

When approaching a cascading tipping point, the long-range cross-correlation between the leading state x and the following state y is expected to increase. The state x becomes more auto-correlated and is less susceptible to noise, and therefore influences y through the coupling in a more robust way. To find long-range cross-correlations, the so-called “detrended cross-correlation analysis” (DCCA) method was developed (Podnobik and Stanley, 2008; Zebende, 2011; Zhang et al., 2001; Zhou, 2008). Instead of using the auto-covariance to calculate the fluctuation function, as is used in Eq. (8), one uses the cross-covariance as follows:

$$\begin{array}{ll}{\displaystyle }{F}_{\mathrm{DCCA}}^{\mathrm{2}}(\mathit{\nu },s)={\displaystyle \frac{\mathrm{1}}{s}}\sum _{i=\mathrm{1}}^s& {\displaystyle }\left[\right(X\left(\right(\mathit{\nu }-\mathrm{1})s+i)-x_{\mathit{\nu }}\left(i\right))\\ \text{(11)}& {\displaystyle }& {\displaystyle }\cdot \left(Y\right((\mathit{\nu }-\mathrm{1})s+i)-y_{\mathit{\nu }}(i\left)\right){]}^{\mathrm{2}}.\end{array}$$

The symbols in Eq. (11) are similar to Eq. (8). With this function, one can calculate the fluctuation function and subsequent power-law scaling coefficient (Podnobik and Stanley, 2008; Zhang et al., 2001), similar to Eq. (9).

A variation on this method was proposed by Zebende (2011) and involves the ratio between $F_{\mathrm{DCCA}}^{\mathrm{2}}$ and F_DFA of the two systems. Specifically, one chooses a certain segment size s and computes

$$\begin{array}{}\text{(12)}& {\displaystyle }{\mathit{\rho }}_{\mathrm{DCCA}}={\displaystyle \frac{F_{\mathrm{DCCA}}^{\mathrm{2}}}{F_{\mathrm{DFA}\mathit{\left\{}x\mathit{\right\}}}{F}_{\mathrm{DFA}\mathit{\left\{}y\mathit{\right\}}}}},\end{array}$$

which measures the level of the long-term cross-correlation between variable x and y; the quantity ρ_DCCA has values between −1 and 1.

There are several a priori limitations of using the power-law scaling coefficient and ρ_DCCA. First of all, the results are sensitive to choices of the maximum segment size (s) and the window size, induced by the noise in our simulations. Making the window size too large decreases the possibility of seeing any development prior to the tipping points, as windows including the tipping event itself are biased by strong auto-correlation and the (tipping) trend in the data. However, making the window size too small increases the uncertainty in the exponential fit. Similarly, adding small segments co-determining the exponential fit makes the method prone to noise, while larger segments are limited by the window size and computation time. In our simulations, window sizes of 120 data points, and segments sizes between 10 and 60 within those windows were used to calculate the fluctuations per segment size and to calculate F(s) (Eq. 9). More research is needed to find the optimal values for the window and segment sizes, and to potentially reduce this limitation in the analysis of cascading transitions. Another limitation of using the power-law scaling coefficient and ρ_DCCA is that the way that the two systems are coupled (e.g. linearity, with or without temporal lag) affects how an evolution in the leading system affects its cross-correlation with the following system in both magnitude and time. Finally, to observe trends in these metrics, the signal in the cross-correlation between the two systems has to overcome the (partly mutually independent) noise. However, close to critical transition, the recovery from noise actually decreases, making the often subtle “change” in the detrended cross-correlation harder to distinguish from noise. These limitations may make the detrended cross-correlation metrics less useful in applications, but trends in the detrended cross-correlation metrics might still act as early warnings for cascading transitions, as shown in the next section.

3.3 Simulations

In this section we discuss the previously described metrics applied in ensemble simulations of cascading transition events. This provides insight into the statistical characteristics of these events, the causal relation between tipping of both systems and the potential prediction of these events. We focus on the double-fold and the fold–Hopf cascading tipping cases for multiple reasons. First, these cases are most illustrative in terms of relation to physical systems. Second, in these cases the leading system starts with a fold bifurcation, which creates a clear threshold for the start of the event (which is convenient for analysis purposes).

Figure 3 Ensemble (100-member) simulations of a dynamical system undergoing a double-fold cascade (Eq. 13) where both systems undergo a transition, parameter values as in Table 2; (a) states of x (black) and y (red), (b) variance of x (black) and y (red), (c) auto-regressive coefficient at lag 1 of x (black) and y (red), (d) detrended fluctuation analysis scaling exponent of x (black) and y (red), (e) detrended cross-correlation analysis scaling exponent and (f) detrended cross-correlation coefficient by Zebende (2011). White-shaded areas indicate windows containing the actual transitions. The increased variance, AR(1) and the DFA scaling exponent prior to transition in the respective leading system and following system, confirms the predicted increased memory through critical slowing down.

Table 2 Parameter values, coupling and initial conditions for the ensemble simulations of the double Fold and fold–Hopf systems as shown in Figs. 3, 4 and 5.

Figure 4 As in Fig. 3, but without any transition in the following system. In this case, only the leading system has a transition. Parameter values are given in Table 2.

3.3.1 Double-fold cascading tipping

To simulate these events and use statistical indicators, noise has to be included. The system of equations used here is

$$\begin{array}{}\text{(13)}& {\displaystyle }\left\{\begin{array}{l}\frac{\mathrm{d}x}{\mathrm{d}t}=a_{\mathrm{1}}{x}^{\mathrm{3}}+a_{\mathrm{2}}x+\mathit{\varphi }+{\mathit{\zeta }}_x\\ \frac{\mathrm{d}y}{\mathrm{d}t}=b_{\mathrm{1}}{y}^{\mathrm{3}}+b_{\mathrm{2}}y+\mathit{\gamma }\left(x\right)+{\mathit{\zeta }}_y,\end{array}\right.\end{array}$$

where in addition to Eq. (3), ζ_x,ζ_y are Gaussian white noise terms. We simulate an ensemble of 100 members with the parameter settings and initial conditions as displayed in Table 2. The results of this ensemble are displayed in Fig. 3. Running windows containing the transition are shaded white as these data are misleading when one wants to know what happens before the bifurcation points. We make the distinction between the leading-transitional period (LTP), which is the time series before the tipping point in the leading system, and the following-transitional period (FTP), which is the time series between the first tipping point and the tipping point in the following system. The FTP can be seen as an in-between state that is stable for both systems, although the following system is close to bifurcation, meaning it is rather prone to noise. Therefore, the duration of this state is highly unpredictable. However, as we will see in this section, the FTP does provide information regarding the diagnosis of a potential second transition.

In the LTP, we clearly see the gradual increasing leading system's variance, the AR(1) coefficient and DFA scaling coefficient. These are all evidence of the leading system slowly approaching a bifurcation point, according to the change in the parameter. There is not much evidence of long-range auto-correlations in the time series of the following system, as its variance is low and the DFA scaling exponent remains below 0.5, pointing towards the fact that the detrended fluctuations are statistically white noise. The AR(1) coefficient of the following system does increase just prior to the first tipping, but it also stays small (compared to unity).

The detrended cross-correlation scaling exponent (abbreviated here as DXA) does give >0.5 values, but the range throughout the ensemble members is too large most of the time to see any structural development when approaching the bifurcation point. During the leading system's transition, a strong increase is visible, pointing towards rather strong cross-correlated behaviour during this period (as the following system also shifts its equilibrium a little).

The quantity ρ_DCCA seems to attain a small positive value (around 0.3) and stays relatively constant throughout the whole time interval. One important aspect of the calculation of ρ_DCCA, as we found by experimentation, is that the values are very sensitive to the segment size s and the moving window size. The moving window determines the amount of data that is available to find long-range correlations, and the segment size has a strong impact on the accuracy of the fits and on the related segmented fluctuations. Due to the fact that we need a temporal evolution of the statistical indicators in this type of problem, we need moving windows and, thus, encounter a problem. As these indicators (DXA and ρ_DCCA) have been successfully applied in simpler systems (Podnobik and Stanley, 2008; Zebende, 2011; Zhang et al., 2001; Zhou, 2008), more research on the sensitivity of the indicators with respect to the segment size and moving window size may lead to more robust results.

During the FTP, the variance, AR(1) and DFA of the leading system are strongly reduced. However, the gradual increase of the following system's variance, AR(1) coefficient and DFA scaling coefficient are definitely visible, pointing towards the approach of a bifurcation in the following system. Also notable is the contrast in the DFA of the following system before and after the tipping of x. The DFA of y goes from a white-noise regime (around 0.5) before the tipping of x towards a regime where the detrended fluctuations point to long-range auto-correlations after the tipping of x (1–1.5). This illustrates the relation between the leading system's state and the following system's DFA scaling exponent. The DXA remains relatively high, but overall no structural development is seen in this graph. The quantity ρ_DCCA exhibits the same behaviour as in the LTP, most likely due to the previously mentioned reasons.

Table 3 Comparison of the ratios of auto-regressive variables prior to and after the first transition, using the ensembles shown in Figs. 3 and 4.

To assess the consequences of the cascading effect on the statistics mentioned, we compare the results with a case where the system (Eq. 13) does not undergo a tipping in the following system (so only one transition remains). The resulting ensemble results are shown in Fig. 4. The most important differences between a regular cascading event and a single tipping event can be found when comparing the variance, AR(1) and DFA scaling coefficient changes between LTP and FTP (or the period after the first transition). During the LTP, the leading system is close to transition and has strong auto-regressive behaviour, whilst the following system is far from its bifurcation point. During the FTP, the following system generally gains memory because it is brought closer to its transition point, and the leading system loses memory because it just arrived at a new state. Therefore, we expect that from the LTP towards the FTP, the variance, AR(1) and DFA will decrease in the leading system, and increase in the following system.

Table 4 Results of Student's t test on the differences between the ratios (in Table 3) of the cases with and without second tipping.

^*p-values calculated using the scipy.stats Python package.

To quantify this effect, Table 3 shows the ratios of the different quantities, indicated by Q, during the FTP and LTP phases (${\stackrel{\mathrm{‾}}{Q}}_{\mathrm{FTP}}/{\stackrel{\mathrm{‾}}{Q}}_{\mathrm{LTP}}$), for the cases “with” a second tipping (corresponding to runs shown in Fig. 3) and “without” second tipping (Fig. 4). All ensemble members are included in these values, accounting for a mean and standard deviation of the ratios. As expected, the leading system's auto-regressive metrics decrease in both cases, which is visible in the fact that the mean values of the ratios of the leading system's auto-regressive variables are lower than 1. Furthermore, as expected, the following system's auto-regressive behaviour increases (ratios >1) in both cases; however, it is striking that in the case of a cascading tipping event (with second tipping), the following system's ratios are much higher than those in the case of a single tipping event (without second tipping). To investigate whether the difference in the ratios between single or double tipping is indeed significant, a Student's t test was carried out. The results are shown in Table 4. The high p values for the leading system's ratios indicate no significant difference between single or double tipping, but the low p values for the following system's ratios indicate a significant difference. This shows the potential of using the ratio of auto-regressive variables before and after a transition to assess whether a cascading transition may follow. Further research is needed to quantify this expectation and to assess the sensitivity of these ratios to the system's parameters.

Figure 5 As in Fig. 3, but for the fold–Hopf cascade (Eq. 14). Parameter values are given in Table 2.

3.3.2 Fold–Hopf cascading tipping

Many statistical indicators have been specifically applied on fold bifurcations, because these transitions show a clear sign of critical slowing down and increased auto-correlation; this is due to their irreversibility and the process of going from one equilibrium towards another. A state approaching a Hopf bifurcation reacts differently to perturbations than a state approaching a saddle-node bifurcation. Therefore, we will consider the fold–Hopf cascade in the light of the previously described statistical indicators. For this, we use the following stochastic dynamical system:

$$\begin{array}{}\text{(14)}& {\displaystyle }\left\{\begin{array}{l}\frac{\mathrm{d}x}{\mathrm{d}t}=a_{\mathrm{1}}{x}^{\mathrm{3}}+a_{\mathrm{2}}x+\mathit{\varphi }+{\mathit{\zeta }}_x\\ \frac{\mathrm{d}y}{\mathrm{d}t}=b_{\mathrm{1}}z+b_{\mathrm{2}}\left(\mathit{\kappa }\right(x)-(y^{\mathrm{2}}+z^{\mathrm{2}}\left)\right)y+{\mathit{\zeta }}_y\\ \frac{\mathrm{d}z}{\mathrm{d}t}=c_{\mathrm{1}}y+c_{\mathrm{2}}\left(\mathit{\kappa }\right(x)-(y^{\mathrm{2}}+z^{\mathrm{2}}\left)\right)z+{\mathit{\zeta }}_z.\end{array}\right.\end{array}$$

This system is similar to Eq. (4), but white noise is added through the terms ζ_x, ζ_y and ζ_z. We used an ensemble of 100 simulations with the parameter settings and initial conditions as displayed in Table 2. The results of the ensemble are shown in Fig. 5. Here, we do not make the distinction between the LTP and the FTP, because in contrast to the double-fold cascade, the following system undergoes a transition that is easily reversed and the system is either stationary or oscillating. We subtract a running average from the states and calculate the statistics from those series to prevent the oscillation from dominating the signal in the auto-correlation. The following system (red) is quickly drawn towards the equilibrium state $(y,z)=(\mathrm{0},\mathrm{0})$, and the leading system (black) is in a steady state. During the time leading up to the bifurcation point, the variance, AR(1) coefficient and DFA of the leading system x gradually increase, as is expected as we force this system towards its bifurcation point. It also seems that the AR(1) of the following system slightly increases during this period.

After the transition of the leading system, the oscillation of the following system immediately starts due to noise. The variance and AR(1) of the following system after the bifurcation are strongly increased with respect to before the bifurcation, despite the removal of the running average to eliminate the oscillation's signal. On average, the DFA scaling exponent of the following system also increases, which relates it to the leading system's state. The DXA sharply increases just prior to the critical transition, although it retains relatively high values throughout the whole time series. The reason behind this might be the low level of noise that is chosen, or other simulation-specific parameters. It could also due to the fact the following system has a high, weakly varying DFA scaling exponent on average, which in turn might affect the height and variability of the cross-correlation. The ρ_DCCA coefficient remains positive and small, as in the double-fold case. Again, this may have to do with the choice of window and segment sizes.

4 Application: the coupled AMOC–ENSO system

In this section, we present an application of the concept of cascading tipping (the fold–Hopf case). This application reflects the fact that cascading transitions are not only a purely mathematical concept, but also occur in idealized physical models. Here, we consider cascading tipping in a model that couples the Atlantic meridional overturning circulation (AMOC) and the El Niño–Southern Oscillation (ENSO).

4.1 Background

To demonstrate and quantify the coupling between the AMOC and ENSO, we use output from global climate model simulations. In the ESSENCE project (Ensemble SimulationS of Extreme weather events under Nonlinear Climate changE) several simulations were performed using the ECHAM5/MPI-OM coupled climate model, including so-called “hosing” experiments Sterl et al. (2008), where fresh water is added around Greenland to mimic ice sheet melting.

Figure 6 (a) Evolution of the five standard SRES-A1b runs (black) and five HOSING-1 runs (red) in terms of the overturning stream function. (b) NINO3.4 of the standard SRES-A1b ensemble (black) and the HOSING-1 ensemble (red). Shaded thin lines indicate monthly means, thick lines indicate the deseasonalised values.

Table 5 NINO3.4 statistics (of deseasonalised data) for the different ensembles. The uncertainty stated is the standard deviation among the five runs within the ensemble. It is visible that in the case of a collapsed overturning, El Niño intensifies more than without a collapsed overturning. Values in bold represent increases in the variability of NINO3.4.

From these climate model simulations we used two ensembles; the first is the “standard” experiment, where greenhouse gases evolve according to observations and from the year 2000 onwards follow the SRES-A1b scenario (experiment name SRES-A1b). The second ensemble is the same as the standard experiment but has additional freshwater input (1 Sv =10⁶ m³ s⁻¹) around Greenland from the end of the year 2000 onwards (experiment name HOSING-1). The HOSING-1 ensemble contains a classical hosing experiment, following the procedure from Jungclaus et al. (2006). Five runs of each ensemble are chosen, specifically runs 041–045 of the HOSING-1 and runs 021–025 of the SRES-A1b ensemble. The temporal resolution used is monthly data between 1950 and 2100. The spatial fields are on a curvilinear grid, with 40 vertical levels in the ocean. We use deseasonalised data because we are interested in interannual variability, not in seasonal variability, as El Niño is associated with these timescales. We use the maximum of the Atlantic meridional overturning stream function at 35^∘ N as an AMOC index, and the NINO3.4 index as an ENSO index, which is the average SST over the region 170–120^∘ W by 5^∘ S–5^∘ N.

The results for the evolution of the AMOC and ENSO are shown in Fig. 6. It is clearly visible that the AMOC decreases strongly in the hosing experiments, by approximately 85 %. Table 5 compares statistical properties for the time intervals before and after 2001 (which is the year at which the hosing starts). We use the non-anomaly statistics, as this gives us information about the differences in the mean. We do note that we only use five runs per ensemble, which means that the uncertainty is not statistically robust. We only present it in Table 5 to give the reader an idea of the range of the variables among the different runs.

It is visible from Table 5 that the variability of NINO3.4 increases (bold numbers) if we compare the 1950–2000 and 2001–2100 periods. This increased variability is visible in both the standard and the HOSING-1 runs. However, the variability is much stronger in the HOSING-1 experiment, indicating that the decrease of the AMOC indeed has an amplifying effect on ENSO. The large difference between the standard and hosing runs suggests that the NINO3.4 index changed in the hosing experiment, potentially as a consequence of the decrease of the AMOC. As the first tipping is artificially induced (without any measurable critical slowing down), and the models used here are far more complex than the simple dynamical systems we discussed in previous sections, the question of whether cascading tipping is actually occurring in these runs is beyond the scope of this paper. We only use this data to justify the coupling of AMOC and ENSO.

Several mechanisms have been suggested in the literature regarding the coupling between the AMOC and ENSO. The first mechanism is concerned with oceanic waves. A colder North Atlantic creates density anomalies that induce the southward propagation of oceanic Kelvin waves (along the American coast) across the Equator. In western Africa, this energy radiates as Rossby waves towards the north and south, which induces Kelvin waves to move along the tip of southern Africa into the Indian Ocean, which eventually reach the Pacific. Consequently, the eastern equatorial Pacific thermocline deepens on a decadal timescale. This deepening has a weakening effect on the amplitude of ENSO (Timmermann et al., 2005).

The second mechanism is concerned with the trade winds. Cooling in the northern tropical Atlantic (due to AMOC weakening) induces anticyclonic atmospheric circulation (Xie et al., 2007) that intensifies the northerly trade winds over the northeastern tropical Pacific. This leads to a southward displacement of the Pacific intertropical convergence zone (Zhang and Delworth, 2005) and generates a meridional SST anomaly due to anomalous heat transport and the wind–evaporation–SST feedback in the Pacific. Also, Dong and Sutton (2007) found an atmospheric coupling through Rossby waves sent into the northeast tropical Pacific. This is in line with Dijkstra and Neelin (1995), who argued that part of the contribution to the zonal wind stress, τ_ext, arises from processes outside the tropical Pacific. The result of the wind-stress coupling between the two systems is an intensification of ENSO, and this mechanism is argued to be stronger than the coupling through oceanic waves (Timmermann et al., 2005).

4.2 Models and coupling

To study the possible cascading transition through the wind-stress coupling, we use a conceptual model. For the AMOC, the classical Stommel box model (Stommel, 1961) is used. It consists of a polar (subscript p) and an equatorial box (subscript e), both with a temperature T and salinity S coupled by a density-driven flow rate. The state variables are then defined as $\mathrm{\Delta }T=T_{\mathrm{e}}-T_{\mathrm{p}}$ and $\mathrm{\Delta }S=S_{\mathrm{e}}-S_{\mathrm{p}}$. The time evolution of these variables is as follows (Cessi, 1994):

$$\begin{array}{}\text{(15)}& {\displaystyle }\left\{\begin{array}{l}\frac{\mathrm{d}\mathrm{\Delta }T}{\mathrm{d}t}=-\frac{\mathrm{1}}{t_{\mathrm{r}}}(\mathrm{\Delta }T-{\mathit{\theta }}_{\mathrm{0}})-Q\left(\mathrm{\Delta }\mathit{\rho }\right)\mathrm{\Delta }T\\ \frac{\mathrm{d}\mathrm{\Delta }S}{\mathrm{d}t}=\frac{F_s}{H}{S}_{\mathrm{0}}-Q\left(\mathrm{\Delta }\mathit{\rho }\right)\mathrm{\Delta }S,\end{array}\right.\end{array}$$

where the first term in the temperature equation refers to relaxation towards a background temperature, and the second term refers to density-driven meridional transport. Specifically, t_r is the surface temperature restoring timescale and θ₀ is the equator-to-pole atmospheric temperature difference. Q(Δρ) is the transport function, which is calculated from a diffusion timescale and the meridional density gradient Δρ. In the salinity equation, S₀ is a reference salinity, and H is the ocean depth. The parameter F_s is the freshwater flux, which can be used as a bifurcation parameter. The stream function

$$\mathrm{\Psi }={\mathit{\gamma }}_{\mathrm{0}}\mathrm{\Delta }\mathit{\rho }/{\mathit{\rho }}_{\mathrm{0}}={\mathit{\gamma }}_{\mathrm{0}}({\mathit{\alpha }}_T\mathrm{\Delta }T-{\mathit{\alpha }}_s\mathrm{\Delta }S)$$

represents the strength of the AMOC, with γ₀>0 a flow constant, ρ₀ a reference density and α_T, α_S the thermal and haline expansion/contraction coefficients.

For the El Niño–Southern Oscillation, we use the conceptual model as proposed in Timmermann et al. (2003). This model has a state vector consisting of the temperature of the western Pacific T₁, the temperature of the eastern Pacific T₂ and the thermocline depth of the western Pacific h₁. The model finds is based on the Zebiak and Cane (1987) ENSO model, with a two-strip and two-box approximation, and a shallow-water model for the upper ocean with a fixed mixed layer depth.

$$\begin{array}{}\text{(16)}& {\displaystyle }\left\{\begin{array}{l}\frac{\mathrm{d}{T}_{\mathrm{1}}}{\mathrm{d}t}=-\mathit{\alpha }(T_{\mathrm{1}}-T_{\mathrm{r}})-\frac{u(T_{\mathrm{2}}-T_{\mathrm{1}})}{L/\mathrm{2}}\\ \frac{\mathrm{d}{T}_{\mathrm{2}}}{\mathrm{d}t}=-\mathit{\alpha }(T_{\mathrm{2}}-T_{\mathrm{r}})-\frac{w(T_{\mathrm{2}}-T_{\mathrm{sub}})}{H_{\mathrm{m}}},\end{array}\right.\end{array}$$

where 1∕α is a typical thermal damping timescale, T_sub is the temperature below the mixed layer, H_m and L are the depths of the mixed layer and the basin width, respectively, w is upwelling velocity and u is atmospheric zonal surface wind linear to wind stress: $u/(L/\mathrm{2})=\mathit{ϵ}\mathit{\beta }\mathit{\tau }$ and $w/H_{\mathrm{m}}=-\mathit{\zeta }\mathit{\beta }\mathit{\tau }$. The parameters ϵ and ζ refer to the strength of zonal and vertical advection (bifurcation parameters).

The subsurface temperature T_sub is parameterized as

$$\begin{array}{}\text{(17)}& T_{\mathrm{sub}}=T_{\mathrm{r}}-{\displaystyle \frac{T_{\mathrm{r}}-T_{\mathrm{r}\mathrm{0}}}{\mathrm{2}}}[\mathrm{1}-{\displaystyle \frac{\tanh (H+h_{\mathrm{2}}-z_{\mathrm{0}})}{h^{*}}}],\end{array}$$

where h₂ is the east equatorial Pacific thermocline depth (calculated as the deviation from a reference depth H), z₀ is the depth for which w becomes its characteristic value and h^* is the sharpness of the thermocline. The thermocline depths are calculated as follows:

$$\begin{array}{}\text{(18)}& {\displaystyle }\left\{\begin{array}{l}{h}_{\mathrm{2}}=h_{\mathrm{1}}+bL\mathit{\tau }\\ \frac{\mathrm{d}{h}_{\mathrm{1}}}{\mathrm{d}t}=r(-h_{\mathrm{1}}-\frac{bL\mathit{\tau }}{\mathrm{2}}),\end{array}\right.\end{array}$$

where b is the efficiency of the wind stress τ to drive the thermocline tilt. For further details, we refer to Timmermann et al. (2003). In the Stommel–Timmermann models, we use the standard parameter settings, as given in the references, unless stated otherwise.

Figure 7 Zonal equatorial wind stress versus the Atlantic temperature gradient. Data from the ESSENCE (Ensemble SimulationS of Extreme weather events under Nonlinear Climate changE) project are used, where the black dots refer to five members of the standard ensemble with SRES-A1b forcing (1950–2100 period ), and the red dots refer to five members of the HOSING-1 ensemble where in 2000 a freshwater perturbation is applied (i.e. 2000–2100 period). A 5-year running mean is applied, yearly averages are shown. The zonal equatorial wind stress is defined here as the average zonal wind stress over the 0–10^∘ N latitudinal band. The Atlantic temperature gradient is defined as the difference between the SST in a northern box (50–60^∘ N, 50–20^∘ W) and a southern box (0–20^∘ N, 45–20^∘ W). The blue line indicates a linear fit.

We couple the AMOC and ENSO models through the relation between the Atlantic meridional temperature gradient (in the Stommel model) and the Pacific zonal wind stress (in the Timmermann model); the Pacific zonal wind stress mechanism, which is the more important of the two abovementioned mechanisms, is found in the literature, and described in the previous section. Even in a simplified model, the relation between wind stress and meridional temperature gradient is physically justified: thermal wind balance indicates a direct connection between the adjustment of wind stress to changes in the meridional temperature gradient. In the Timmermann model, the zonal wind stress τ is expressed as

$$\begin{array}{}\text{(19)}& \mathit{\tau }={\displaystyle \frac{\mathit{\mu }(T_{\mathrm{2}}-T_{\mathrm{1}})}{\mathit{\beta }}},\end{array}$$

with the μ∕β parameters that control the influence of the zonal temperature gradient on the wind stress set to 0.02 Pa⋅ K⁻¹. We add an external wind stress term τ_ext that is dependent on the meridional temperature gradient in the Atlantic ΔT, i.e.

$$\begin{array}{}\text{(20)}& {\displaystyle }\mathit{\tau }={\mathit{\tau }}_{\mathrm{ext}}\left(\mathrm{\Delta }T\right)+{\displaystyle \frac{\mathit{\mu }}{\mathit{\beta }}}(T_{\mathrm{2}}-T_{\mathrm{r}}),\end{array}$$

with a negative relation between τ_ext and the Atlantic meridional SST gradient ΔT as we know from the literature described above (a stronger positive ΔT results in stronger easterlies and thus a negative τ_ext). Note that both the Pacific wind stress τ and specifically its external part τ_ext should always be negative. The total wind stress is negative because this area (at low altitude) is strictly dominated by easterly winds, and τ_ext is negative because, through the meridional temperature gradient, it reflects the influence of the zonal mean Hadley cell on the equatorial Pacific. Physically, the Hadley cell only induces negative zonal wind stress in this region.

In the coupling (Eq. 21), we fix β and vary μ as the coupling parameter. For τ_ext we take a linear relation as follows:

$$\begin{array}{}\text{(21)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{ext}}={\mathit{\alpha }}_{\mathit{\tau }}\mathrm{\Delta }T+{\mathit{\gamma }}_{\mathit{\tau }}-{\mathit{\tau }}_{\mathrm{0}},\end{array}$$

where all coefficients are constant over time. The parameters α_τ and γ_τ can be estimated from the ESSENCE data as discussed in Sect. 4.1; τ₀ reflects the constant part in the zonal mean wind stress, which we subtract because we are interested in the contribution of changes in the meridional overturning. Using five ESSENCE runs per ensemble for both the standard forcing and hosing-experiment, respectively, ΔT is computed as the absolute difference between the SST in the North Atlantic area (50–60^∘ N × 50–20^∘ W) and the equatorial Atlantic region (0–20^∘ N × 45–20^∘ W). For the wind stress τ_ext, the zonally integrated wind stress averaged over the 0–10^∘ N region is taken. In Fig. 7, 5-year running means of annual averages are plotted for the hosing simulations (in red) and the standard simulations (in black). Clearly, τ_ext decreases as ΔT increases, meaning that when the AMOC collapses (larger ΔT) the wind stress τ_ext becomes more negative and the external part of the trade winds increases. However, we also note that the spread in the simulation data is large, which can, in part, be attributed to internal variability present in the simulations. The coefficients α_τ and γ_τ were found to be (from a linear fit) −0.000376 Pa ⋅ ^∘C⁻¹ and −0.0119 Pa. By looking at the ΔT regime in Fig. 7, τ₀ is chosen to be the wind stress at 19 ^∘C: ${\mathit{\tau }}_{\mathrm{0}}={\mathit{\alpha }}_{\mathit{\tau }}\cdot \mathrm{19}+{\mathit{\gamma }}_{\mathit{\tau }}\approx -\mathrm{0.0190}$ Pa. This results in a final quantized expression for the coupling:

$$\begin{array}{}\text{(22)}& {\mathit{\tau }}_{\mathrm{ext}}\approx -\mathrm{0.000376}\cdot \mathrm{\Delta }T+\mathrm{0.00715}.\end{array}$$

Figure 8 Bifurcation diagrams and forward runs of the Stommel (a, b) and Timmermann (c, d) models. The blue arrow indicates the collapse of the overturning circulation, which (negatively) amplifies the external zonal wind stress in the Pacific τ_ext, such that the system enters an oscillatory state. The orange arrow indicates subsequent tipping in the following (ENSO) system. (a) Meridional temperature gradient equilibria versus freshwater flux, and (b) non-dimensional stream function versus freshwater flux. These figures show the multiple states of the overturning. (c) Eastern equatorial Pacific SST versus τ_ext (for μ=0.00145), showing a regime where the system is stationary and a regime where the system is oscillatory; (d) eastern equatorial Pacific SST versus μ for different values of τ_ext. Orange dots indicate Hopf bifurcation points, and orange dotted lines indicate oscillatory regimes. Black and grey solid lines indicate stable equilibria, and black dashed lines indicate unstable equilibria.

4.3 Results

The AMOC model's bifurcation diagrams are shown in Fig. 8a and b, and clearly display a saddle-node structure. For an interval of values of the freshwater flux F_s the system has multiple equilibria, whilst for other values only one equilibrium remains. This means that when we are in the high-Ψ branch and F_s is large enough, the system can make a transition to the low-Ψ branch. This is depicted by the blue arrow in Fig. 8b.

Figure 9 Simulation run of the coupled Stommel–Timmermann model for different model configurations, where the collapse of the overturning flow function (black) leads to the crossing of a Hopf bifurcation in the eastern equatorial Pacific SST (orange). Parameter values as in Timmermann et al. (2003), with μ=0.00146.

The bifurcation diagram of the ENSO model with τ_ext as a parameter is shown in Fig. 8c. The bifurcation diagrams become much simpler than in the original Timmermann et al. (2003) model, the reason for this is extensively discussed in Dijkstra and Neelin (1995). Figure 8d shows the influence of μ for the position of the oscillatory regime: on each branch, two Hopf bifurcations can be found and the μ value of the first Hopf bifurcation decreases with more negative τ_ext. This indicates that a Hopf bifurcation can be crossed if τ_ext is decreased, while μ is kept constant. In other words, for the right value of μ, the eastern Pacific SST starts oscillating (El Niño “intensifies”) when the easterly external wind is increased. For the coupled model, we use μ=0.00146.

Using τ_ext to couple the AMOC and ENSO models, we performed simulations with Δt=0.1 days and the Runge–Kutta fourth-order integration method. To initiate the collapse of the overturning, a freshwater forcing F_s is applied in the form of a step function:

$$\begin{array}{}\text{(23)}& {\displaystyle }{F}_s=\left\{\begin{array}{l}\mathrm{0.006}\text{if}t\le \mathrm{500}y\\ \mathrm{0.01}\text{if}t>\mathrm{500}y.\end{array}\right.\end{array}$$

Using the coupling of Eq. (23), we attain the results shown in Fig. 9. The exact quantification of the coupling partly modulates which effect the collapse of the AMOC has on ENSO. For the chosen coupling, the collapse of the overturning leads to the crossing of the first Hopf bifurcation point in the following system, and an oscillation starts growing. As is visible in Fig. 6, the relation between ΔT and τ has quite some spread, which implies a large uncertainty in the values of α_τ and γ_τ. We would like to stress that the regime in which the Hopf bifurcation is crossed is dependent on multiple variables, such as these coupling parameters. Running the forward integration of the coupled model for values between ${\mathit{\alpha }}_{\mathit{\tau }}=-\mathrm{0.00041}$ and ${\mathit{\alpha }}_{\mathit{\tau }}=-\mathrm{0.00033}$, uncovers the fact that (ceteris paribus), for higher (lower) values of α_τ, the oscillation indeed becomes weaker (stronger), finally resulting in a disappearance of the oscillation at ${\mathit{\alpha }}_{\mathit{\tau }}\approx -\mathrm{0.000335}$ at a time step of 0.25 days.

Despite the parameter sensitivity, this is a typical illustration of the fold–Hopf cascading behaviour discussed in earlier sections. This re-enforces the possibility that cascading transitions are possible in real physical systems.

5 Summary, discussion and conclusions

In this paper, we introduced the concept of cascading tipping, which can occur when a transition in a leading system alters background conditions for a following system that then undergoes a transition. We presented a mathematical framework around this concept, using generic bifurcations (saddle-node and Hopf) in both leading and following systems. Four types of deterministic dynamical systems with the possibility for cascading events were formulated, including the double-fold cascade, the fold–Hopf cascade, the Hopf–fold cascade and the double-Hopf cascade. In all cases we assumed a linear coupling between the following and leading system. The double-fold coupled system has previously been considered in another context (Brummitt et al., 2015), where it was also noted that not all subsystems undergo tipping (“hopping”) in systems with more than two coupled fold cascades. Moreover, stochastically coupled multi-stable systems have been considered in networks, where different types of domino effects can occur depending on the synchrony of the transition in the different network nodes (Ashwin et al., 2017; Creaser et al., 2018). Here we only consider two coupled systems, but allow different types of bifurcations, and the systems are physically coupled in a directional way.

We discussed statistical indicators and analysis tools for cascading tipping points. Indicators for cascading tipping points are found in detrended cross-correlation analysis (DCCA) and a special case of extrapolation using the DFA of the following system. These tools were applied in simulations involving both the double-fold and fold–Hopf cascades. The increased variance, AR(1) and DFA scaling exponent are clearly found in each case of single tipping. The cross-correlation indicators (DCCA and ρ_DCCA) did not evolve much throughout the time series, which indicates their insensitivity with respect to the proximity to single tipping points. Several limitations on the use of these variables have been mentioned. However, it seems that these metrics are highly sensitive to window and segment sizes, which means that their potential as early warnings of cascading transition events is inconclusive. The ratios of auto-regressive metrics before and after the first transition seem to be a stronger warning of cascading transitions. More research is needed to exactly quantify these metrics.

The concept of cascading tipping was applied to study the behaviour of a model describing a link between the Atlantic meridional overturning circulation (MOC) and ENSO. We modelled this using a coupling between the Stommel (1961) model and the Timmermann et al. (2003) ENSO model by a meridional temperature gradient-dependent term in the external wind stress of the ENSO model. Through analysis of the bifurcation diagrams and simulations, a cascading tipping event is indeed possible within this model in the form of the fold–Hopf cascade. Obviously, both models are highly idealized and more detailed models of both AMOC and ENSO are needed to demonstrate the occurrence of such a cascading transition in the climate system.

A potential example of a double-fold cascade, which was not further treated here, could be the impact of a bistable MOC on the (bistable) land ice formation on the Antarctic continent. In this case the coupling exists through the atmospheric CO₂ concentration, which depends on mixing and circulation in the ocean while strongly determining the existence of an ice sheet DeConto and Pollard (2003). During the Eocene–Oligocene transition, where a large ice sheet grew on Antarctica, a two-step signal is observed in the deep-sea δ¹⁸O ratio, suggesting two abrupt transitions. Using a box model by Gildor and Tziperman (2000), Tigchelaar et al. (2011) showed that a two-step signal can be produced by (first) a MOC transition which changes the CO₂ concentration meaning that a transition occurs in the land-ice model. Although from a physical perspective, this is a potential example of a cascading transition, we make no claim about whether such a transition likely occurred at the Eocene–Oligocene transition. Here, more detailed models are also needed and transition are expected to be more complicated (Tantet et al., 2018).

These two applications reflect the relevance of this paper. There are likely many cases in which these cascading events occur in climate and therefore highlight the importance of the topic. Future research will point out whether these events are likely to happen in the future climate and whether these effects also occur in fields other than climate science. Of course, this paper covers the very basics of deterministic cascading events. However, one can imagine a wide range of phenomena if more complicated transitions between attractors are considered and when noise is included. For example, when a leading chaotic system is coupled to a deterministic following system with a saddle-node bifurcation structure, a slight change in the chaotic attractor may change the background conditions for the following system such that it undergoes a transition. An application here may be the effect of a mid-latitude atmospheric jet on the Atlantic MOC. We hope that this paper will stimulate more research on the various types of cascading tipping and also on the development of well-suited indicators and early warnings of such events.