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 ﬁnd that there can be a signiﬁcant 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


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).
Published by Copernicus Publications on behalf of the European Geosciences Union.
J. Lohmann et al.: Abrupt climate change as a rate-dependent cascading tipping point 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 iceocean 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., 2009Scheffer et al., , 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.

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 where the density is given by the equation of state of seawater with some reference densities, temperatures and salinities represented by ρ 0 , T 0 and S 0 , respectively.The two model variables represent the dimensionless temperature difference T = α T (T e − T p ) and salinity difference S = α S (S e − S p ) in between the boxes.This defines the dimensionless overturning circulation strength Temperature and salinity in the boxes relax towards an atmospheric temperature and salinity forcing T a e,p and S a e,p .The meridional difference of the forcing drives the circulation and is represented by the two parameters η 1 ∝ (T a e − T a p ) and η 2 ∝ (S a e − S a p ).A third parameter represents the timescale ratio of the temperature and salinity relaxation η 3 = τ T τ S .The model is then defined by the stochastic differential equations 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 η 1 and η 2 .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 nonsmooth 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 bistability depends on the parameter η 3 .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.

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 a p , thus modulating the temperature forcing gradient η 1 ∝ (T a e − T a p ).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  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 = −0.282.We define the state at R = 0 as the stadial state, yielding a fixed point with positive sea ice cover I + 0 ≈ 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 icecovered 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).
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 η 1 , which we now define as η 1 (I ) = η 0 1 −κ (I )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: 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 η 0 1 = 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 η 1 = 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 η 1 .Here we choose a scenario where during the stadial the sea ice reduces the atmospheric forcing from η 1 = 3.0 to η 1 = 2.65.κ is then chosen such that η 1 = 2.65 at the stadial fixed point I + 0 and η 1 = 3.0 for I < 0, yielding κ = 0.35/I + 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 η 1 (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.

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 (σ T = σ S = 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 η 1 and η 2 as control parameter, it occurs for a larger range of amplitudes and rates of the parameter shift when changing η 1 .
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 η 1 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 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 dyhttps://doi.org/10.5194/esd-12-819-2021 Earth Syst.Dynam., 12, 819-835, 2021 namics 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.
The critical ramping duration depends on the amplitude η 1 of the parameter shift (Fig. 7a).Rate-induced tipping becomes possible at a certain η 1 , where the basin instability condition is first satisfied.Increasing η 1 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.

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 η 1 = 3.0 for higher noise levels.
By introducing noise, tipping becomes a mixture of rateinduced 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.
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 dehttps://doi.org/10.5194/esd-12-819-2021 Earth Syst.Dynam., 12, 819-835, 2021 terministic 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 σ T = σ S = 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.

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 η 1 (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 η 1 (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 η 1 (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 η 1 (I ) from η 1 = 2.65 (at the stadial sea ice fixed point for R = 0) to η 1 = 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 shows trajectories for a cascading stadialinterstadial transition in the deterministic limit when R is ramped down from R = 0 to R = −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 Earth Syst.Dynam., 12, 819-835, 2021 https://doi.org/10.5194/esd-12-819-2021tries 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 η 1 (I ) towards η 1 = 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 oceanonly 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 https://doi.org/10.5194/esd-12-819-2021 Earth Syst.Dynam., 12, 819-835, 2021 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.

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, 2011or 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.

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 uncer- tainty 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.
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 https://doi.org/10.5194/esd-12-819-2021 Earth Syst.Dynam., 12, 819-835, 2021 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 J , as defined in Eq. (A6). Figure 15a-f shows that J 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 J .The figure also shows probabilities p of observing a value of J estimated around the attractor that is larger than a value of J in the vicinity of the saddle.This measures the performance of J 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.
An operational early warning signal can be constructed by estimating J in a sliding window, and raising an alert as soon as a threshold J c is exceeded.Choosing a location of J 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 J comes close to a perfect classifier, detecting the saddle approach with very low false positive and very high true positive rates.Figure 16 shows J estimated from time series in a sliding window, along with critical slowing down indicators.J 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.

Discussion
In this work we propose a conceptual model that describes a mechanism for abrupt climate change comprising a rateinduced 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 https://doi.org/10.5194/esd-12-819-2021 Earth Syst.Dynam., 12, 819-835, 2021 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 rateinduced 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.

Conclusions
Building on previous studies of proxy records and state-ofthe-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.

Figure 2 .
Figure 2. Schematic of the coupled sea ice-ocean model including model parameters and variables (bold).A description of the parameters is given in Table1.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 ∝ T e − T p and S ∝ S e − S p .In the polar ocean box, the sea ice cover I insulates the ocean from the cold atmospheric temperature T a p .

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

Figure 4 .
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 η 1 as control parameter with η 2 = 1.0 and η 3 = 0.3.(c, d) The same as (a, b) but with η 2 as control parameter with η 1 = 3.0 and η 3 = 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).

Figure 5 .
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 η 1 = 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 η 1 = 3.0.Two trajectories where η 1 is ramped linearly from η 1 = 2.65 to η 1 = 3.0 within 300 and 500 years are shown in blue and purple, respectively.The initial conditions T , S = (2.4,2.5) are indicated by the yellow cross.Open symbols indicate the positions of the fixed points at η 1 = 2.65, and the black curve indicates the corresponding basin boundary from (a).

Figure 6 .
Figure 6.Time series of q = T − S in the Stommel model when ramping the parameter from η 1 = 2.65 to η 1 = 3.0 at different rates.The realizations are initialized at T , S = (2.4,2.5), which is close to the off fixed point at η 1 = 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 η 1 = 3.0 are shown as horizontal lines.

Figure 7 .
Figure 7. (a) Critical ramping duration below which there is a rate-induced tipping in the Stommel model when shifting the parameter from η 1 = 2.65 to η 1 = 2.65 + η 1 .The off attractor loses stability in the bifurcation at η off 1 = 3.333, as indicated by the dashed red line.(b) Normalized shortest distance to the basin boundary B as a function of the normalized distance to the bifurcation λ = (η off 1 − η 1 ) • (η off 1 − η on 1 ) −1 .η on 1 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 .
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 η 1 = 2.65 to η 1 = 3.00.Different noise levels σ T = σ S = σ 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.

Figure 9 .
Figure 9. (a-c) Three realizations in phase space of the Stommel model with σ T = σ S = 0.2, where η 1 is ramped from η 1 = 2.65 to η 1 = 3.00 over 300 years.The filled dot (triangle) marks the off fixed point (saddle) at η 1 = 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 = (2.4,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 = 0.1 used to define the time of tipping is shown as a dotted line.

Figure 10 .
Figure 10.Cascading stadial-interstadial transition in the coupled sea ice-ocean model where R is ramped from R = 0 to R = −0.3within 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 η 1 (I (t)) = η 0 1 − κ (I (t)) • I (t).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 η 1 (I (t)).

Figure 12 .
Figure 12.Ensemble simulations of the coupled sea ice ocean model, where R is ramped linearly from R = 0 to R = −0.3within 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. (cf) 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 .
Figure 13.(a) Simulation in phase space of the Stommel model with σ T = σ S = 0.2, where η 1 is ramped from η 1 = 2.65 to η 1 = 3.00 within 300 years.The two dotted lines correspond to the levels q = T − S = 0.06 and q = T − S = 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 .

Figure 14 .
Figure 14.Distributions of variance and autocorrelation for ensembles of time series from the Stommel model (σ T = σ S = 0.2).These are estimated from time series segments around the initial fixed point at η 1 = 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 η 1 = 3.0 (dotted).

Figure 15 .
Figure 15.(a-f) Distributions of the early warning indicator J for ensembles of time series from the Stommel model (σ T = σ S = 0.2), estimated around the initial fixed point at η 1 = 2.65 (black) and close to the saddle point (orange).For each realization, J 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 J 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.

Figure 16 .
Figure 16.Early warning indicators estimated in a 200-year sliding window from an ensemble of time series of the Stommel model, where η 1 is ramped from η 1 = 2.65 to η 1 = 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 J (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.

Table 1 .
Description of model parameters.