the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Statistical precursor signals for Dansgaard–Oeschger cooling transitions
Niklas Boers
Given growing concerns about climate tipping points and their risks, it is important to investigate the capability of identifying robust precursor signals for the associated transitions. In general, the variance and shortlag autocorrelations of the fluctuations increase in a stochastically forced system approaching a critical or bifurcationinduced transition, making them theoretically suitable indicators to warn of such transitions. Paleoclimate records provide useful test beds if such a warning of a forthcoming transition could work in practice. The Dansgaard–Oeschger (DO) events are characterized by millennialscale abrupt climate changes during the glacial period, manifesting most clearly as abrupt temperature shifts in the North Atlantic region. Some previous studies have found such statistical precursor signals for the DO warming transitions. On the other hand, statistical precursor signals for the abrupt DO cooling transitions have not been identified. Analyzing Greenland ice core records, we find robust and statistically significant precursor signals of DO cooling transitions in most of the interstadials longer than roughly 1500 years but not in the shorter interstadials. The origin of the statistical precursor signals is mainly related to socalled rebound events, humps in the temperature observed at the end of interstadial, some decades to centuries prior to the actual transition. We discuss several dynamical mechanisms that give rise to such rebound events and statistical precursor signals.
 Article
(4819 KB)  Fulltext XML

Supplement
(1797 KB)  BibTeX
 EndNote
A tipping point is a critical threshold beyond which a system reorganizes, often abruptly and/or irreversibly (IPCC, 2023). Once a tipping point is passed, a system can abruptly transition to an alternative stable or oscillatory state (Boers et al., 2022). Empirical and modeling evidence suggests that some components of the Earth system might indeed exhibit tipping behavior, which poses arguably one of the greatest potential risks in the context of ongoing anthropogenic global warming (Armstrong McKay et al., 2022; Boers et al., 2022). Paleoclimate evidence supports the fact that abrupt climate changes due to crossing tipping points actually occurred in the past (Dakos et al., 2008; Brovkin et al., 2021; Boers et al., 2022). The Dansgaard–Oeschger events (Dansgaard et al., 1993) are one of such past abrupt climate changes during the last glacial period and the focus of this study.
Tipping point behavior is mathematically classified into three different types (Ashwin et al., 2012). (1) Bifurcationinduced tipping is an abrupt or qualitative change of a system owing to a bifurcation of a stable state (more generally a quasistatic attractor). (2) Noiseinduced tipping is an escape from a neighborhood of a quasistatic attractor caused by the action of noisy fluctuations (Ditlevsen and Johnsen, 2010). (3) Rateinduced tipping occurs when a system fails to track a continuously changing quasistatic attractor (Ashwin et al., 2012; Wieczorek et al., 2023; O'Sullivan et al., 2023). In realworld systems, tipping behaviors often result from a combination of several of the above (Ashwin et al., 2012).
The theory of critical slowing down (CSD) provides a framework to anticipate critical (or bifurcationinduced) transitions (Carpenter and Brock, 2006; Scheffer et al., 2009; Kuehn, 2013; Boers, 2018, 2021; Boers and Rypdal, 2021; Boers et al., 2022; Bury et al., 2020). The framework is based on the fact that the stability of a stable state is gradually lost as the system approaches the bifurcation point. Theoretically, the variance of the fluctuations around the fixed point diverges and the autocorrelation with a sufficiently small lag increases toward 1 at the critical point of a codimension1 bifurcation (Boers et al., 2022; Scheffer et al., 2009; Bury et al., 2020), where the codimension1 bifurcations are, in simple terms, the bifurcations that can be typically encountered by the change of a single control parameter (Thompson and Sieber, 2011). Thus, the changes in CSD indicators such as the increase of the variance as well as the autocorrelation can be seen as statistical precursor signals (SPSs) of critical transitions. See Dakos et al. (2012) as well as Boers et al. (2022) for various methods and CSD indicators for anticipating critical transitions.
Dansgaard–Oeschger (DO) events are millennialscale abrupt climate transitions during glacial intervals (Dansgaard et al., 1993). They are most clearly imprinted in the δ^{18}O and calcium ion concentration [Ca^{2+}] records from the Greenland ice cores (Fig. 1) (Rasmussen et al., 2014; Seierstad et al., 2014). The δ^{18}O and [Ca^{2+}] are interpreted as proxies for site temperature and atmospheric circulation changes, respectively. DO warmings occur typically within a few decades and are followed by gradual cooling during relatively warm glacial states termed “interstadials”, before a rapid return to cold states referred to as “stadials”. The amplitude of the abrupt warming transitions ranges from 5 to 16.5 °C (Kindler et al., 2014, and references therein). The Greenland temperatures change concurrently with the North Atlantic temperatures (Bond et al., 1993; Martrat et al., 2004), atmospheric circulation patterns (Yiou et al., 1997), seawater salinity (Dokken et al., 2013), seaice cover (Sadatzki et al., 2019), and the Atlantic Meridional Overturning Circulation (AMOC), as inferred from indices such as $\mathrm{Pa}/\mathrm{Th}$ and δ^{13}C (Henry et al., 2016). The combined proxy evidence suggests that the DO events arise from interactions among these components (Menviel et al., 2020; Boers et al., 2018). The prevailing view is that the mode switching of the AMOC plays a principal role in generating DO events (Broecker et al., 1985; Rahmstorf, 2002), but it remains debated whether the inferred AMOC changes are a driver of DO events or a response to the changes in the atmosphere–ocean–seaice system in the North Atlantic, Nordic Seas, and the Arctic (Li and Born, 2019; Dokken et al., 2013). Recently, an increasing number of comprehensive climate models succeeded in simulating DOlike selfsustained oscillations, suggesting that DO events can arise spontaneously from complex interactions between the AMOC, ocean stratification/convection, atmosphere, and sea ice (Peltier and Vettoretti, 2014; Vettoretti et al., 2022; Brown and Galbraith, 2016; Klockmann et al., 2020; Zhang et al., 2021; Kuniyoshi et al., 2022; MalmiercaVallet and Sime, 2023).
The DO events are considered the archetype of climate tipping behavior (Boers et al., 2022; Brovkin et al., 2021). Early works found an SPS based on autocorrelation for one specific DO warming – the onset of Bølling–Allerød interstadial (Dakos et al., 2008). In subsequent works, the existence of SPS for DO warmings was questioned considering that DO warmings are noise induced rather than bifurcation induced (Ditlevsen and Johnsen, 2010; Lenton et al., 2012). However, a couple of later studies detected SPS for several DO warmings either by ensemble averaging of CSD indicators for many events (Cimatoribus et al., 2013) or by using wavelettransform techniques focusing on a specific frequency band (Rypdal, 2016; Boers, 2018). On the other hand, it has so far not been shown whether DO coolings are preceded by characteristic CSDbased precursor signals as well.
Recent studies have inferred that the AMOC is currently at its weakest in at least a millennium (Rahmstorf et al., 2015; Caesar et al., 2018) (see also Kilbourne et al. (2022) for possible uncertainties). The declining AMOC trend is projected to continue in the coming century, although the projections of the AMOC strength in the next 100 years are model dependent (MassonDelmotte et al., 2021). Furthermore, the studies applying CSD indicators to observed AMOC fingerprints (Boers, 2021; BenYami et al., 2023; Ditlevsen and Ditlevsen, 2023) as well as a longterm reconstruction of the Atlantic multidecadal variability (Michel et al., 2022) suggest that the AMOC stability may have declined and the AMOC might be approaching a dangerous tipping point. In this context, it is important to investigate whether CSDbased precursor signals can be detected for the DO cooling transitions as well, supposing that the DO events reflect past AMOC changes. However, of course, predictability of past events does not necessarily imply any predictability in the future, especially given that the recent AMOC weakening is presumably driven by global warming and is thus from a mechanistic point of view different from past AMOC weakenings in the glacial period.
In this study we report SPS for DO cooling transitions recorded in δ^{18}O and log _{10}[Ca^{2+}] (Seierstad et al., 2014; Rasmussen et al., 2014) from three Greenland ice cores: the North Greenland Ice Core Project (NGRIP), the Greenland Ice Core Project (GRIP), and the Greenland Ice Sheet Project 2 (GISP2) (see Fig. 1 for NGRIP). The important source of observed SPS stems from socalled rebound events, humps in the temperature proxy occurring at the end of interstadials, some decades to centuries prior to the transition (Capron et al., 2010). When CSD indicators such as variance or lag1 autocorrelation are used for anticipating a tipping point, we conventionally assume that a system gradually approaches a bifurcation point. However, if DO cycles are spontaneous oscillations as suggested in the comprehensive climate models (see above), in a strict sense there might not be any bifurcation occurring around the timings of the abrupt transition in the DO cycles. Nevertheless, with conceptual models of DO cycles, we demonstrate that CSD indicators may show precursor signals for abrupt transitions due to particular dynamics near a critical point or a critical manifold.
The remainder of this paper is organized as follows. In Sect. 2, the data and method are described. In Sect. 3, we identify robust and statistically significant SPSs for several DO cooling transitions following interstadials with sufficient data length and show that rebound events prior to cooling transition are the source of observed SPSs. In Sect. 4 we discuss the results by using conceptual models. A summary is given in Sect. 5.
2.1 Greenland ice core records
We explore CSDbased precursor signals for DO cooling transitions recorded in δ^{18}O and log _{10}[Ca^{2+}] (Seierstad et al., 2014; Rasmussen et al., 2014) from three Greenland ice cores: NGRIP, GRIP, and GISP2 (see Fig. 1 for NGRIP). Multiple records are used for a robust assessment because each has regional fluctuations as well as proxy and icecoredependent uncertainties. The six records have been synchronized and are given at 20year resolution (Seierstad et al., 2014; Rasmussen et al., 2014). They continuously span the last 104 kyr b2k (kiloyears before 2000 CE), beyond which only NGRIP δ^{18}O is available up to a part of the Eemian Interglacial. In addition, we use a version of the NGRIP δ^{18}O and dust records at 5 cm depth resolution (EPICA community members, 2004; Gkinis et al., 2014; Ruth et al., 2003) in order to check the dependence of results on temporal resolutions, with the caveat that these highresolution records span only the last 60 kyr.
We follow the classification of interstadials and stadials and associated timings of DO warming and cooling transitions by Rasmussen et al. (2014), where Greenland interstadials (stadials) are labeled with “GI” (“GS”) with a few exceptions below. A rebound event is an abrupt warming often observed before an interstadial abruptly ends (Capron et al., 2010) (arrows in Figs. 1, 2, and 3). Generally, a long interstadial accompanies a long rebound event (their durations are correlated with R^{2}=0.95, Capron et al., 2010). In Capron et al. (2010) and Rasmussen et al. (2014), GI14 and the subsequent interstadial, GI13, are seen as one long interstadial, with GI13 considered to be a strongly expressed rebound event ending GI14 because the changes in δ^{18}O and log _{10}[Ca^{2+}] during the quasistadial GS14 do not reach the baseline levels of adjacent stadials. Similarly GI23.1 and GI22 are also seen as one long interstadial, and GI22 is regarded as a rebound event (and GS23.1 as quasistadial) (Capron et al., 2010; Rasmussen et al., 2014). GI20a is also recognized as a rebound event in Rasmussen et al. (2014). Given that the rebound events are warmings following a colder spell during interstadial conditions that does not reach the stadial levels (Rasmussen et al., 2014), we regard the following nine epochs as reboundtype events: GI8a, the hump at the end of GI11 (42 240–∼42 500 yr b2k), GI12a, GI13, the hump at the end of GI16.1 (56 500–∼56 900 yr b2k), GI20a, GI21.1cba (two warming transitions), GI22, and GI25a. When we examine the effect of rebound events on our results, we exclude the entire parts including the cold spells prior to the rebound events.
The start (warming) and end (cooling) of each DO event are identified in 20year resolution based on both δ^{18}O and [Ca^{2+}] in Rasmussen et al. (2014). The estimated uncertainty of event timing varies from event to event. We remove the 2σ uncertainty range of the event timing (40 to 400 years) estimated in Rasmussen et al. (2014) from our calculation of CSD indicators. It effectively excludes parts of the transitions themselves from the calculation of the CSD indicators. Since the calculation of CSD indicators requires a minimum length of data points, we mainly focus on interstadials longer than 1000 years after removing 2σ uncertainty ranges of the transition timings, using 20year resolution data (Table S1 in the Supplement). This results in 12 DO interstadials to be investigated (Fig. 1, gray shades). We deal with the interstadials shorter than 1000 years but longer than 300 years using highresolution records in Sect. 3.2.
2.2 Statistical indicators of critical slowing down
Based on the theory of critical slowing down (CSD), we posit that the stability of a dynamical system perturbed by noise is gradually lost as the system approaches a bifurcation point (Boers et al., 2022; Scheffer et al., 2009; Bury et al., 2020). For the fold bifurcation (also known as the saddlenode bifurcation), the variance of the fluctuations around a local stable state diverges and the autocorrelation function of the fluctuations increases toward 1 at any lag τ. The same is true for the transcritical as well as the pitchfork bifurcation (Bury et al., 2020). For the Hopf bifurcation, the variance increases, but the autocorrelation function of the form $C\left(\mathit{\tau}\right)={e}^{\mathit{\nu}\left\mathit{\tau}\right}\mathrm{cos}\mathit{\omega}\mathit{\tau}$ may increase or decrease depending on τ, where ν(≤0) and ±ωi are the real and imaginary parts of the complex eigenvalues of the Jacobian matrix of the local linearized system (Bury et al., 2020). Nevertheless, the autocorrelation function C(τ) increases for sufficiently small τ. For discrete time series, we follow previous studies and calculate the lag1 autocorrelation corresponding to a minimal τ. These characteristics can be used to anticipate abrupt transitions caused by codimension1 bifurcations. Promisingly previous studies show that these CSD indicators and related measures can indeed anticipate simulated AMOC collapses (Boulton et al., 2014; Klus et al., 2018; Livina and Lenton, 2007; Held and Kleinen, 2004).
Prior to calculating CSD indicators, we estimate the local stable state by using a local regression method called the locally weighted scatterplot smoothing (LOESS) (Cleveland et al., 1992; Dakos et al., 2012). In this approach, the time series is seen as a scatter plot and fitted locally by a polynomial function, giving more weight to points near the point whose response is being estimated and less weight to points further away. Here the polynomial degree is set to 1; i.e., the smoothing is performed with the local linear fit. Nevertheless, the LOESS provides nonlinear smoothed curves. The smoothing span parameter α that defines the fraction of data points involved in the local regression is set to 50 % of each interstadial length in a demonstration case, but we examine the dependence of results on α over the range 20 %–70 %. The difference between the record and the smoothed one gives the residual fluctuations. The CSD indicators, i.e., variance and lag1 autocorrelation, are calculated for the residuals over a rolling window. The size W of this rolling window is set to 50 % of each interstadial length in the demonstration case. To test the robustness, this is changed over the range 20 %–60 %.
The statistical significance of precursor signals of critical transitions, in terms of positive trends of CSD indicators, is assessed by hypothesis testing (Theiler et al., 1992; Dakos et al., 2012; Rypdal, 2016; Boers, 2018). We consider as null model a stationary stochastic process with preserved variance and autocorrelation. The n surrogate data are prepared from the original residual series by the phaserandomization method, thus preserving the variance and autocorrelation function of the original time series via the Wiener–Khinchin theorem. Here we take n=1000. The linear trend (a_{o}) of the CSD indicator for the original time series and the linear trends (a_{s}) of CSD indicators for the surrogate data are calculated. We consider the precursor signal of the original time series statistically significant at the 5 % level if the probability of a_{s}>a_{o} (p value) is less than 0.05. The significance level of 0.05 is conservative given that some works analyzing ecological or paleodata adopt the significance level of 0.1 (Dakos et al., 2012; Thomas et al., 2015).
3.1 Characteristic precursor signals of DO coolings
As CSD indicators we consider the variance and lag1 autocorrelation, calculated in rolling windows across each interstadial. The 12 interstadials longer than 1000 years are magnified in Figs. 2 and 3 (top rows, blue) for the NGRIP δ^{18}O record. See Figs. S1–S10 in the Supplement for the other records. For each interstadial, the nonlinear trend is estimated using LOESS smoothing (Figs. 2 and 3, top row, red). In this case the smoothing span α that defines the fraction of data points involved in the local regression is set to 50 % of each interstadial length. Gaussian kernel smoothing gives similar results. The difference between the record and the nonlinear trend gives the approximately stationary residual fluctuations (second row). The CSD indicators are calculated from the residual series over a rolling window. In Figs. 2 and 3 the rolling window size W is set to 50 % of each interstadial length (a default value in Dakos et al., 2008). The smoothing span α and the rolling window size W are taken as fractions of individual interstadial length because timescales of local fluctuations (such as the duration of rebound events) change with the entire duration of interstadial. We examine the dependence of the results on α and W as part of our robustness tests.
The variance is plotted in the third row of Figs. 2 and 3. Positive trends in the variance are observed for 9 out of 12 interstadials; the individual trends are statistically significant in 6 out of 12 cases (p<0.05), based on a null model assuming the same overall variance and autocorrelation, constructed by producing surrogates with randomized Fourier phases. The lag1 autocorrelation is also plotted for the same data in the bottom row. Positive trends in lag1 autocorrelation are observed for 10 out of 12 interstadials but are statistically significant only in two cases (p<0.05). Just a positive trend without significance cannot be considered a reliable SPS, but if one indicator has a significantly positive trend, the other indicator with a consistently positive trend may at least support the conclusion (e.g., GI19.2 and GI1413 in Fig. 3). In several cases (GI24.2, GI21.1, GI16.1, GI1413, and GI12), the lag1 autocorrelation first decreases and then increases. The initial decreases, harming monotonic increases of CSD indicators, might reflect a memory of the preceding DO warming transition. On the other hand, the drastic increases in both indicators near the end of the interstadials reflect the rebound events (arrows in Figs. 2 and 3). We obtain similar results for the other ice core records (Figs. S1–S10). While we observe a number of positive trends for all the records, the number of statistically significant trends detected depends on the record and CSD indicator (Table S2).
We check the robustness of our results against changing smoothing span α and rolling window size W (Dakos et al., 2012). We calculate the p value for the trend of each indicator changing the smoothing span between 20 % and 70 % of interstadial length (in steps of 10 %) and the rolling window size between 20 % and 60 % (also in 10 % steps). This yields a 6×5 matrix for the p values. Example results for GI25 and δ^{18}O are shown in Fig. 4a (variance) and 4b (lag1 autocorrelation). The cross mark (x) indicates significant positive trends (p<0.05) and the small open circle (o) indicates positive trends that are significant at the 10 % level but not at the 5 % level. Full results for the 12 interstadials, six records, and two CSD indicators are shown in Figs. S11–S22. We consider positive trends in CSD indicators, i.e., the SPS of the transition, to be overall robust if we obtain significant positive trends (p<0.05) for more than half (>15) of the 30 parameter sets.
The robustness analysis is performed for all the long interstadials of the six records and the two CSD indicators (Fig. 4d). Among the 12 interstadials, we find at least one robust SPS for eight interstadials (GI25, GI23.1, GI21.1, GI20, GI19.2, GI14, GI12, and GI8) and multiple robust SPSs for six (GI25, GI23.1, GI21.1, GI14, GI12, and GI8). If the data series is a stationary stochastic process, the probability of spuriously observing a robust SPS is estimated to be 5 % (Appendix A). In this case, the probability of detecting more than two robust SPSs from 12 independent stationary samples (i.e., from each row in Fig. 4d) by chance is only ∼2 %. Thus, the risk of obtaining the results by chance is quite low. For each interstadial, the detection of robust SPSs depends on the proxy and core. This is possibly because different proxies from different cores are contaminated by different types and magnitudes of noise (e.g., δ^{18}O may record local fluctuations of temperatures and log _{10}[Ca^{2+}] turbulent fluctuations of local wind circulations). Robust SPSs are observed for most interstadials longer than roughly 1500 years (GI25, GI23.1, GI21.1, GI20, GI19.2, GI14, GI12, and GI8 except GI1) but not for the other interstadials, shorter than roughly 1500 years (compare Figs. 4c and 4d).
3.2 Further sensitivity analyses
We examine how much the rebound events affect the detection of CSDbased SPS. For this purpose CSD indicators are again calculated excluding the rebound events and their preceding cold spells (see Sect. 2.1). While eight interstadials (GI25, GI23.1, GI21.1, GI20, GI16, GI14, GI12, and GI8) exhibit robust SPS with the rebound events included, only four interstadials (GI23.1, GI14, GI12, and GI8) exhibit robust SPS without the rebound events (Fig. S23). The rebound events should hence be considered important, sometimes indispensable, sources for SPS of DO coolings.
We also examine the dependence of the results on the time resolution of the data. Here we use a highresolution (5 cm depth) δ^{18}O record (EPICA community members, 2004; Gkinis et al., 2014) and a dust record (Ruth et al., 2003) from the NGRIP over the last 60 kyr. Since the data in these records are nonuniform in time, they are linearly resampled every 5 years before calculating CSD indicators. We focus on 11 interstadials longer than 300 years in order to have enough data points. For the dust record, three interstadials (GI15, GI8, and GI7) are excluded from the analysis because the original data have long parts of missing values. The CSD indicators, calculated with a smoothing span of α=50 % and rolling windows of W=50 %, are shown in Figs. S24–S27. Through the robustness analyses with respect to α and W, we find at least one robust SPS for three out of 11 interstadials (Fig. S28). The robust SPSs for GI1413 and GI12 from the highresolution records are consistent with those from the 20year resolution records. Moreover for GI1, the highresolution δ^{18}O record exhibits a robust SPS in terms of lag1 autocorrelation, although the 20year resolution record does not. Robust SPSs have not been observed again for interstadials shorter than roughly 1500 years (Figs. S28 and 4).
We detected robust precursor signals of DO cooling transitions for most interstadials longer than roughly 1500 years but not for shorter interstadials. The results suggest that long interstadials, the existence of rebound events, and the presence of SPS for the DO cooling transitions are all related (except for GI19.2, which has no noticeable rebound event). These aspects may be related to generic properties of nonlinear dynamical systems. On the basis of conceptual mathematical models, we discuss four possible dynamical mechanisms leading to the precursor signals of DO cooling transitions. In three of four mechanisms, oscillations such as the rebound events can systematically arise prior to the abrupt cooling transitions. These modeling results justify the inclusion of the rebound events in the search for precursor signals presented above. Unless otherwise mentioned, details on model parameters as well as the hysteresis experiments conducted are given in Appendix B.

The fold bifurcation mechanism. Since the pioneering work by Stommel (1961), the AMOC is considered to exhibit bistability depending on the background condition (Rahmstorf, 2002). The bistability of the AMOC strength x may be conceptually modeled by a doublefold bifurcation model: $\dot{x}=f\left(x\right)+\mathit{\mu}+\mathit{\xi}\left(t\right)$, where f(x) has two fold points such as x−x^{3} and $\leftx\right(\mathrm{1}x)$. Here we take the quadratic from $f\left(x\right)=\leftx\right(\mathrm{1}x)$, but the following arguments are qualitatively the same for x−x^{3}. The parameter μ represents the surface salinity flux (i.e., negative freshwater flux), and ξ(t) denotes white Gaussian noise. The unperturbed model for ξ(t)=0 has equilibria on an Sshaped curve: $f\left(x\right)+\mathit{\mu}=\mathrm{0}$ (Fig. 5a, green). The state x(t) initially on the upper stable branch jumps down to the lower stable branch as μ decreases across the fold bifurcation point at $\mathit{\mu}=\mathrm{0.25}$. The variance and the autocorrelation of the local fluctuations (i.e., CSD indicators) increase as μ slowly approaches the fold bifurcation point since the restoring rate toward the stable state decreases, as shown in Fig. 6a (Scheffer et al., 2009; Boers et al., 2022).

Stochastic slow–fast oscillation mechanism. The FitzHugh–Nagumo (FHN) system is a prototypical model for slow–fast oscillations and excitability (FitzHugh, 1961; Nagumo et al., 1962). It is often invoked for conceptual models of DO oscillations (Rial and Yang, 2007; Kwasniok, 2013; Roberts and Saha, 2017; Mitsui and Crucifix, 2017; Lohmann and Ditlevsen, 2019; Riechers et al., 2022; Vettoretti et al., 2022). An FHNtype model of DO oscillations can be obtained by introducing a slow variable y into the fold bifurcation model: $\dot{x}=\frac{\mathrm{1}}{{\mathit{\tau}}_{x}}\left(\rightx\left\right(\mathrm{1}x)+y+\mathit{\mu})+\mathit{\xi}\left(t\right)$, $\dot{y}=\frac{\mathrm{1}}{{\mathit{\tau}}_{y}}(xy)$, where τ_{x} and τ_{y} are timescale parameters (τ_{x}≪τ_{y}). Invoking the saltoscillator hypothesis for DO oscillations suggested by the comprehensive climate model simulations that are successful in reproducing DO cycles (Vettoretti and Peltier, 2018), we may interpret y as the surface mixedlayer salinity in the northern North Atlantic and Labrador Sea, which gradually decreases (increases) when the AMOC intensity x is strong (weak).
Here we consider the case that the system is excitable. For example, for μ=0.26, the unperturbed system has a stable equilibrium near the upper fold point of the Sshaped critical manifold, $\mathit{\left\{}\right(x,y)\in {\mathbb{R}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}y=\leftx\right(\mathrm{1}x)\mathit{\mu}\mathit{\}}$ (Fig. 5c, green), but the dynamical noise ξ(t) enables the escape from the barely stable equilibrium and sustains stochastic oscillations (Fig. 5b and c, blue). Due to the timescale separation (τ_{x}≪τ_{y}), the oscillations occur along the attracting parts of the critical manifold (Fig. 5c). Because y is much slower than x, the dynamics of x is similar to the dynamics of the fold bifurcation model with slowly changing y. As a result, SPS can be effectively observed near the upper fold point of the critical manifold (Fig. S29). However, this example is not rigorous bifurcationinduced tipping. In the example of an excitable system (Fig. 5b and c), the underlying system always has a weakly stable fixed point, and no true bifurcation leading to CSD occurs. In fact, the actual tipping in this case is noise induced. However, we can effectively observe the SPS in CSD indicators in this case as well, since the system would in each cyclic iteration move from more stable to less stable conditions until it finally tips to initiate the next cycle, and this partial decrease in stability is imprinted in the CSD indicators (Fig. S29). The increase of the variance prior to the transitions in the FHN model is also reported in Meisel and Kuehn (2012). Since the unperturbed system has an equilibrium near the upper fold point, the motion is stagnant near the fold point. This provides favorable conditions for observing SPS. The state jumps from the upper branch of the critical manifold to its lower branch often occur after an upward jump induced by noise. These upward jumps resemble the rebound events prior to DO cooling transitions. The overall phenomenon is the same in the selfsustained oscillation regime of the FHN model, as long as the equilibrium is located near the upper fold point of the critical manifold (μ≃0.25).

Hopf bifurcation mechanism. In contrast to the fold bifurcation, the Hopf bifurcation manifests oscillatory instability (Strogatz, 2018). In several ocean box models, the thermohaline circulation is destabilized via a Hopf bifurcation (Alkhayuon et al., 2019; Lucarini and Stone, 2005; Abshagen and Timmermann, 2004; Sakai and Peltier, 1999). It is also considered a potential generating mechanism of DO oscillations in a loworder coupled climate model (Timmermann et al., 2003) and in a comprehensive climate model (Peltier and Vettoretti, 2014). Assume that the parameter μ decreases slowly in the FHNtype model (Fig. 5d and e). The underlying dynamics changes from the stable equilibrium to the limitcycle oscillations at the Hopf bifurcation point $\mathit{\mu}={\mathit{\mu}}_{\text{Hopf}}\equiv (\mathrm{1}{\mathit{\tau}}_{x}/{\mathit{\tau}}_{y}{)}^{\mathrm{2}}/\mathrm{4}$ (Appendix B). If stochastic forcing is added to the system, noiseinduced small oscillations can appear prior to the onset of the limitcycle oscillations (Fig. 5d and e). The precursor oscillations resemble rebound events, while their shape depends on the noise as well as the change rate of μ. Again SPS can be observed near the Hopf bifurcation point (Fig. S30) (Bury et al., 2020; Meisel and Kuehn, 2012; Boers et al., 2022). The small oscillations prior to downward transitions, like the DO rebound events, do not appear if the system goes deeply into the selfsustained oscillation regime away from the Hopf bifurcation point ($\mathit{\mu}<{\mathit{\mu}}_{\text{Hopf}}\approx \mathrm{0.245}$ in Fig. 5d).

Mixedmode oscillation mechanism. Mixedmode oscillations (MMOs) are periodic oscillations consisting of small and largeamplitude oscillations (Koper, 1995; Desroches et al., 2012; Berglund and Landon, 2012). They often arise in systems with one fast variable and two slow variables (Desroches et al., 2012). In this regard, we can extend the above FHNtype model to exhibit MMOs, for example, as follows: $\dot{x}=\frac{\mathrm{1}}{{\mathit{\tau}}_{x}}\left(\rightx\left\right(\mathrm{1}x)+y+\mathit{\mu})$, $\dot{y}=\frac{\mathrm{1}}{{\mathit{\tau}}_{y}}(xy+k(zy\left)\right)$ and $\dot{z}=\frac{\mathrm{1}}{{\mathit{\tau}}_{z}}(xz+k(yz\left)\right)$, where z is another slow variable with timescale τ_{z} (≫τ_{x}) and k is the diffusivecoupling constant between slow variables. We interpret y as the surface salinity in the northern North Atlantic convection region that directly affects the AMOC strength x again and z as the surface salinity outside the convection region that affects the surface salinity y in the convection region via mixing. This extended FHNtype model is introduced here only to demonstrate that MMOs may appear in an FHNtype model with a minimal dimensional extension. For certain parameter settings (Appendix B), the system has an unstable equilibrium $(x,y,z)=(\sqrt{\mathit{\mu}},\sqrt{\mathit{\mu}},\sqrt{\mathit{\mu}})$ of saddlefocus type, with one stable direction with a negative real eigenvalue and a twodimensional unstable manifold with complex eigenvalues with a positive real part. The slow–fast oscillations occur along the critical manifold $\mathit{\left\{}\right(x,y,z)\in {\mathbb{R}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}y=\leftx\right(\mathrm{1}x)\mathit{\mu}\mathit{\}}$ (Fig. 5f and g). However, due to the saddlefocus equilibrium on the critical manifold, the trajectory is attracted toward the saddle from the direction of the stable manifold (black segment) and repelled from it in a spiralling fashion. The striking point is the systematic occurrence of smallamplitude oscillations prior to the abrupt transition, which also resemble the rebound events prior to the DO cooling transitions. A more realistic time series is obtained if an observation noise is added on x(t) (Fig. S31). Then SPS can be stably observed near the fold point of the critical manifold.
Based on the four types of simple mathematical models, we have proposed four possible dynamical mechanisms for the DO cooling transitions that can manifest statistical precursor signals (SPS): (1) strict CSD due to the approaching of a fold bifurcation; (2) CSD in a wider sense, in stochastic slow–fast oscillations; (3) noiseinduced oscillations prior to Hopf bifurcations; or (4) the signal of mixedmode oscillations. While the details of these mechanisms are different, they are all related to the fold points of the equilibrium curve or the critical manifolds. As a result, the SPS can be detected by the conventional CSD indicators.
Mechanisms (2), (3), and (4) can generate behavior resembling the rebound events, leading to increases in the classic CSD indicators. In the toy models, rebound eventlike behavior is generated when the trajectory passes by an equilibrium point with marginal stability (i.e., the equilibrium has neither strong stability leading to a permanent state nor strong instability leading to short interstadials) (Fig. 5b–g). In this case, the duration of the modeled interstadial is relatively long in relation to the marginal stability. By contrast, the absence of equilibrium or the presence of a strongly unstable equilibrium near the fold point of the critical manifold leads to brief interstadials without a rebound event and consequently a lack of SPS. This provides a possible explanation of why the rebound events and the robust SPS are simultaneously observed for long interstadials but not for short interstadials.
Another possible explanation for the lack of SPS for short interstadials is the following. The common assumption underlying CSD theory is that the parameter change is much slower than the system's relaxation time, and the latter is much slower than the correlation time of the noise (Thompson and Sieber, 2011; Ashwin et al., 2012). If this assumption is violated, it is difficult to detect CSDbased SPS (Clements and Ozgul, 2016; van der Bolt et al., 2021). Consider the foldbifurcationinduced tipping in the Stommeltype model (1), for example (Fig. 6). If the change in the parameter μ is faster than the system's relaxation time toward the moving stable equilibrium, it is unlikely to detect significant CSDbased SPS (Fig. 6b) because the trajectory evolves systematically away from the equilibrium and thus cannot feel the flattening of the potential around the equilibrium even at the true bifurcation point.
In this study we have explored statistical precursor signals (SPSs) and significant increases in critical slowing down (CSD) indicators (variance and lag1 autocorrelation), for Dansgaard–Oeschger (DO) cooling transitions following interstadials, using six Greenland ice core records. Among the 12 interstadials longer than 1000 years, we find at least one robust SPS for eight interstadials longer than roughly 1500 years (GI25, GI23.1, GI21.1, GI20, GI19.2, GI14, GI12, and GI8) and multiple robust SPSs for six of them (GI25, GI23.1, GI21.1, GI14, GI12, and GI8) (Fig. 4d). Robust SPSs are, however, not observed for interstadials shorter than roughly 1500 years. One might link the increase in the proxy variance with the tendency of larger climatic fluctuations in colder climates (Ditlevsen et al., 1996), but the increases in the lag1 autocorrelation cannot generally be explained by it. The analysis removing the rebound events from the data shows that the rebound events prior to the cooling transitions are important in producing the SPS.
We have proposed four different dynamical mechanisms to explain the observed SPSs: (1) strict CSD due to the approaching of a fold bifurcation; (2) CSD in a wider sense, in stochastic slow–fast oscillations; (3) noiseinduced oscillations prior to Hopf bifurcations; or (4) the signal of mixedmode oscillations. In the last three mechanisms, oscillations such as the rebound events can systematically arise prior to the abrupt cooling transitions. These precursor oscillations are due to marginally (un)stable equilibria on the critical manifolds that cause a longlived quasistable state (like long interstadials). This can explain why rebound events and SPSs are simultaneously observed only for long interstadials and are not observed for short ones. While the SPSs for bifurcationinduced tipping events (mechanisms 1 and 3) are established, detailed properties of SPSs for the stochastic slow–fast oscillations of the excitable system (mechanism 2) and for the mixedmode oscillations (mechanism 4) remain to be elucidated.
We should mention the assumptions made in this study as well as alternative scenarios for the DO cooling transitions. First, the four dynamical mechanisms assume slow changes in parameters or slow variables which cause bifurcations in the fast subsystem. On the other hand, the rateinduced tipping mechanism has also been invoked for a possible AMOC collapse, where the rate of change of the external forcing (e.g., freshwater flux or atmospheric CO_{2} concentration) determines the future AMOC state (Alkhayuon et al., 2019; Lohmann and Ditlevsen, 2021; Ritchie et al., 2023). The lack of observed SPSs for the interstadials less than roughly 1500 years indicates a ratedependent aspect of the DO cooling transitions. However, a comprehensive investigation of DO cooling transitions from the viewpoint of rateinduced tipping is beyond the scope of this work. Second, a recent study using an ocean general circulation model shows that a reboundeventtype behavior of AMOC is caused by a behavior called “the intermediate tipping”, due to multiple stable ocean circulation states that exist near but prior to the tipping point leading to a significant AMOC weakening (Lohmann et al., 2023). The intermediate tipping mechanism for rebound events is different from the possible lowdimensional dynamical mechanisms proposed in this study. Further studies are needed to elucidate the dynamical as well as physical origin of DO coolings and associated rebound events.
We have shown that past abrupt DO cooling transitions in the North Atlantic region can be anticipated based on classic CSD indicators if they are preceded by long interstadials. However, it is found to be difficult to anticipate DO cooling events, at least from the 20yearresolution ice core Greenland records, if they occur after a short interstadial. If the DO coolings transitions are actually associated with AMOC weakening (see the “Introduction”), our results may have an implication for the predicted weakening of the AMOC and its possible collapse in the future: the prediction with CSD indicators could be more difficult if the forcing changes fast. There is, however, a caveat to this implication because the past DO cooling transitions are different from the presently inferred AMOC changes. The time resolution (mainly 20 years and additionally 5 years) and the length (mainly >1000 years and additionally >300 years) of the interstadial segment data used in this study are coarser and mostly longer than the annual data used for analyzing AMOC fingerprints during the industrial period (Boers, 2021; BenYami et al., 2023; Ditlevsen and Ditlevsen, 2023) and the last millennium (Michel et al., 2022). More crucially, the revealed predictability of past DO cooling events does not necessarily imply predictability of a potential future AMOC collapse since the recent AMOC weakening, possibly driven by global warming but potentially also part of natural variability, is mechanistically very different from past AMOC weakening in the glacial period.
The statistical significance of precursor signals of critical transitions, in terms of positive trends of CSD indicators, is assessed by hypothesis testing (Theiler et al., 1992; Dakos et al., 2012; Rypdal, 2016; Boers, 2018). We consider a stationary stochastic process with preserved variance and autocorrelation as the null model. The n surrogate data are prepared from the original residual data series by the phaserandomization method, thus preserving the variance and autocorrelation function of the original time series via the Wiener–Khinchin theorem. Here we take n=1000. The linear trend (a_{0}) of the CSD indicator for the original residual time series and the linear trends (a_{s}) of CSD indicators for the surrogate data are calculated. We consider the precursor signal of the original series statistically significant at the 5 % level if the probability of a_{s}>a_{o} (p value) is less than 0.05. Thus, if the original data are already a stationary stochastic process (exhibiting no CSD), one should expect spuriously significant results at a probability of 0.05 by definition. In principle this is independent of the smoothing span α as well as the rolling window size W used for calculating CSD indicators. We consider a precursor signal robust if we find significant cases (p<0.05) for more than half (>15) of 30 combinations of α and W. Then the probability of observing a robust precursor signal can be shown to be 0.05. In order to check this numerically, we generate 5000 surrogates of the original δ^{18}O series of interstadial GI25 and calculate the probability of finding robust precursor signals. The resulting fractions are 0.041 for the variance and 0.039 for the lag1 autocorrelation, which are close to 0.05. For the case of GI12, we obtain 0.038 for the variance and 0.047 for the lag1 autocorrelation, again close to 0.05. These results support the probability of observing a robust precursor signal being 5 % if the data are stationary stochastic processes.
Here we describe specific settings for four conceptual models representing different candidate mechanisms for the DO cooling transitions. Unless otherwise mentioned, the stochastic differential equations below are solved with the Euler–Maruyama method with a step size of 10^{−3}.

The bistability of the AMOC strength x can be conceptually modeled by a doublefold bifurcation model: $\dot{x}=f\left(x\right)+\mathit{\mu}+\mathit{\xi}\left(t\right)$, where f(x) has two fold points; here for f one can use either $f\left(x\right)=x{x}^{\mathrm{3}}$ or $f\left(x\right)=\leftx\right(\mathrm{1}x)$. We take the quadratic function $f\left(x\right)=\leftx\right(\mathrm{1}x)$ that arises in the Stommel (1961) model. μ represents a forcing parameter on the AMOC strength x, e.g., salinity forcing on the North Atlantic (i.e., negative freshwater forcing). ξ(t) is white Gaussian noise, e.g., freshwater perturbations or weather forcing. In Fig. 5a, we set $\sqrt{\langle {\mathit{\xi}}^{\mathrm{2}}\rangle}=\mathrm{0.03}$, and the initial condition is taken at x(0)=1.1, near the upper stable fixed point of the unperturbed system. The parameter μ is then slowly decreased from 0.1 to −0.4 over the period from t=0 to t=500, to trigger the bifurcationinduced transition.

The FitzHugh–Nagumotype (FHNtype) system is a prototypical model of slow–fast oscillators (FitzHugh, 1961; Nagumo et al., 1962) and often invoked for conceptual models of DO oscillations (Rial and Yang, 2007; Kwasniok, 2013; Roberts and Saha, 2017; Mitsui and Crucifix, 2017; Lohmann and Ditlevsen, 2019; Riechers et al., 2022; Vettoretti et al., 2022). The FHNtype model subjected to dynamical noise can be obtained by introducing a slow variable y into the fold bifurcation model: $\dot{x}=\frac{\mathrm{1}}{{\mathit{\tau}}_{x}}\left(\rightx\left\right(\mathrm{1}x)+y+\mathit{\mu})+\mathit{\xi}\left(t\right)$, $\dot{y}=\frac{\mathrm{1}}{{\mathit{\tau}}_{y}}(xy)$, where τ_{x} and τ_{y} are timescale parameters (τ_{x}≪τ_{y}). Following the saltoscillator hypothesis to explain DO cycles (Vettoretti and Peltier, 2018), we may interpret y as the salinity in the polar halocline surface mixed layer, which decreases (increases) when the AMOC is strong (weak). In the case of Fig. 5b and c, we set μ=0.26, τ_{x}=0.01, τ_{y}=1, and $\sqrt{\langle {\mathit{\xi}}^{\mathrm{2}}\rangle}=\mathrm{0.3}$, and the initial condition is taken at $\left(x\right(\mathrm{0}),y(\mathrm{0}\left)\right)=(\mathrm{0.2},\mathrm{0.45})$. The x nullcline (critical manifold) of the unperturbed system is $y=\leftx\right(\mathrm{1}x)\mathit{\mu}$ (Fig. 5c, green) and the y nullcline is $y=x$ (Fig. 5c, magenta dashed). The intersection of the x and y nullclines is the equilibrium point of the unperturbed system $(\sqrt{\mathit{\mu}},\sqrt{\mathit{\mu}})$, which is near the fold point of the critical manifold in this parameter setting.

To demonstrate the Hopf bifurcation mechanism in Fig. 5d and e, the same stochastic FHNtype model is used with τ_{x}=0.01, τ_{y}=1, and $\sqrt{\langle {\mathit{\xi}}^{\mathrm{2}}\rangle}=\mathrm{0.05}$, but here μ is gradually decreased from 0.3 to 0.2, over a period of 5 time units. For $\mathrm{0.2}<\mathit{\mu}<\mathrm{0.3}$, the system has a unique equilibrium point at $(x,y)=(\sqrt{\mathit{\mu}},\sqrt{\mathit{\mu}})$. The Hopf bifurcation of an equilibrium occurs if the complex eigenvalues of the Jacobian matrix at the equilibrium pass the imaginary axis (Strogatz, 2018). The eigenvalues of the Jacobian matrix at this equilibrium are ${\mathit{\lambda}}_{\pm}=\frac{\mathrm{1}}{\mathrm{2}}\mathit{\{}\frac{\mathrm{1}\mathrm{2}\sqrt{\mathit{\mu}}}{{\mathit{\tau}}_{x}}\frac{\mathrm{1}}{{\mathit{\tau}}_{y}}\pm \sqrt{(\frac{\mathrm{1}\mathrm{2}\sqrt{\mathit{\mu}}}{{\mathit{\tau}}_{x}}\frac{\mathrm{1}}{{\mathit{\tau}}_{y}}{)}^{\mathrm{2}}\frac{\mathrm{8}\sqrt{\mathit{\mu}}}{{\mathit{\tau}}_{x}{\mathit{\tau}}_{y}}}\mathit{\}}$. These eigenvalues λ_{±} are complex conjugates for $\frac{\mathrm{1}}{\mathrm{4}}(\mathrm{1}+\frac{{\mathit{\tau}}_{x}}{{\mathit{\tau}}_{y}}\mathrm{2}\sqrt{\frac{{\mathit{\tau}}_{x}}{{\mathit{\tau}}_{y}}}{)}^{\mathrm{2}}<\mathit{\mu}<\frac{\mathrm{1}}{\mathrm{4}}(\mathrm{1}+\frac{{\mathit{\tau}}_{x}}{{\mathit{\tau}}_{y}}+\mathrm{2}\sqrt{\frac{{\mathit{\tau}}_{x}}{{\mathit{\tau}}_{y}}}{)}^{\mathrm{2}}$. In this range of μ, the real part of λ_{±} changes from negative to positive at the Hopf bifurcation point: ${\mathit{\mu}}_{\text{Hopf}}=\frac{\mathrm{1}}{\mathrm{4}}(\mathrm{1}\frac{{\mathit{\tau}}_{x}}{{\mathit{\tau}}_{y}}{)}^{\mathrm{2}}.$ For ${\mathit{\tau}}_{x}/{\mathit{\tau}}_{y}=\mathrm{0.01}$, ${\mathit{\mu}}_{\text{Hopf}}\approx \mathrm{0.245}.$ The initial condition is taken at the origin.

The mixedmode oscillation model is obtained if the FHNtype model is extended to have multiple interacting slow variables. For example, $\dot{x}=\frac{\mathrm{1}}{{\mathit{\tau}}_{x}}\left(\rightx\left\right(\mathrm{1}x)+y+\mathit{\mu}$), $\dot{y}=\frac{\mathrm{1}}{{\mathit{\tau}}_{y}}(xy+k(zy\left)\right)$, and $\dot{z}=\frac{\mathrm{1}}{{\mathit{\tau}}_{z}}(xz+k(yz\left)\right)$, where z is another slow variable with timescale τ_{z} (≫τ_{x}) and k is the diffusive coupling constant between slow variables. We interpret y as the surface salinity in the northern North Atlantic convection region, which directly affects the AMOC strength x, and z as the surface salinity outside the convection region that affects the surface salinity y in the convection region via mixing. We set τ_{x}=0.02, τ_{y}=2, τ_{z}=4, μ=0.225, and k=0.8. This system has an unstable equilibrium $(x,y,z)=(\sqrt{\mathit{\mu}},\sqrt{\mathit{\mu}},\sqrt{\mathit{\mu}})$ of saddlefocus type, with one stable direction with a negative real eigenvalue of −0.67 and a twodimensional unstable manifold with two complex conjugate eigenvalues with a positive real part of 0.94±4.7i. The initial condition is taken at $\left(x\right(\mathrm{0}),y(\mathrm{0}),z(\mathrm{0}\left)\right)=(\mathrm{0.5},\mathrm{0.5},\mathrm{0.5})$.
The Greenland ice core records used in this study can be obtained from https://www.iceandclimate.nbi.ku.dk/data/ (Centre for Ice and Climate, 2024). The R codes used in this study are available from https://doi.org/10.5281/zenodo.10841655 (Mitsui, 2024).
The supplement related to this article is available online at: https://doi.org/10.5194/cp206832024supplement.
TM conceived the study and conducted the analyses with contributions from NB. Both authors discussed and interpreted the results. TM wrote the manuscript with contributions from NB.
The contact author has declared that neither of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
The authors thank Keno Riechers and Maya BenYami for their helpful comments.
The authors acknowledge funding by the Volkswagen Foundation. This is TiPES contribution #243; The TiPES (“Tipping Points in the Earth System”) project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 820970. Niklas Boers acknowledges further funding by the European Union's Horizon 2020 research and innovation programme under the Marie SklodowskaCurie grant agreement no. 956170, as well as from the Federal Ministry of Education and Research under grant no. 01LS2001A.
This paper was edited by Bjørg Risebrobakken and reviewed by three anonymous referees.
Abshagen, J. and Timmermann, A.: An organizing center for thermohaline excitability, J. Phys. Oceanogr., 34, 2756–2760, 2004. a
Alkhayuon, H., Ashwin, P., Jackson, L. C., Quinn, C., and Wood, R. A.: Basin bifurcations, oscillatory instability and rateinduced thresholds for Atlantic meridional overturning circulation in a global oceanic box model, P. Roy. Soc. A, 475, 20190051, https://doi.org/0.1098/rspa.2019.0051, 2019. a, b
Armstrong McKay, D. I., Staal, A., Abrams, J. F., Winkelmann, R., Sakschewski, B., Loriani, S., Fetzer, I., Cornell, S. E., Rockström, J., and Lenton, T. M.: Exceeding 1.5 C global warming could trigger multiple climate tipping points, Science, 377, eabn7950, https://doi.org/10.1126/science.abn7950, 2022. a
Ashwin, P., Wieczorek, S., Vitolo, R., and Cox, P.: Tipping points in open systems: bifurcation, noiseinduced and ratedependent examples in the climate system, Philos. T. Roy. Soc. A, 370, 1166–1184, 2012. a, b, c, d
BenYami, M., Skiba, V., Bathiany, S., and Boers, N.: Uncertainties in critical slowing down indicators of observationbased fingerprints of the Atlantic Overturning Circulation, Nat. Commun., 14, 8344, https://doi.org/10.1038/s41467023440469, 2023. a, b
Berglund, N. and Landon, D.: Mixedmode oscillations and interspike interval statistics in the stochastic FitzHugh–Nagumo model, Nonlinearity, 25, 2303, https://doi.org/10.1088/09517715/25/8/2303, 2012. a
Boers, N.: Earlywarning signals for DansgaardOeschger events in a highresolution ice core record, Nat. Commun., 9, 2556, https://doi.org/10.1038/s41467018048817, 2018. a, b, c, d
Boers, N.: Observationbased earlywarning signals for a collapse of the Atlantic Meridional Overturning Circulation, Nat. Clim. Change, 11, 680–688, 2021. a, b, c
Boers, N. and Rypdal, M.: Critical slowing down suggests that the western Greenland Ice Sheet is close to a tipping point, P. Natl. Acad. Sci. USA, 118, e2024192118, https://doi.org/10.1073/pnas.2024192118, 2021. a
Boers, N., Ghil, M., and Rousseau, D.D.: Ocean circulation, ice shelf, and sea ice interactions explain Dansgaard–Oeschger cycles, P. Natl. Acad. Sci. USA, 115, E11005–E11014, https://doi.org/10.1073/pnas.180257311, 2018. a
Boers, N., Ghil, M., and Stocker, T. F.: Theoretical and paleoclimatic evidence for abrupt transitions in the Earth system, Environ. Res. Lett., 17, 093006, https://doi.org/10.1088/17489326/ac8944, 2022. a, b, c, d, e, f, g, h, i, j
Bond, G., Broecker, W., Johnsen, S., McManus, J., Labeyrie, L., Jouzel, J., and Bonani, G.: Correlations between climate records from North Atlantic sediments and Greenland ice, Nature, 365, 143–147, 1993. a
Boulton, C. A., Allison, L. C., and Lenton, T. M.: Early warning signals of Atlantic Meridional Overturning Circulation collapse in a fully coupled climate model, Nat. Commun., 5, 1–9, 2014. a
Broecker, W. S., Peteet, D. M., and Rind, D.: Does the oceanatmosphere system have more than one stable mode of operation?, Nature, 315, 21–26, 1985. a
Brovkin, V., Brook, E., Williams, J. W., Bathiany, S., Lenton, T. M., Barton, M., DeConto, R. M., Donges, J. F., Ganopolski, A., McManus, J., Praetorius, S., de Vernal, A., AbeOuchi, A., Cheng, H., Claussen, M., Crucifix, M., Gallopín, G., Iglesias, V., Kaufman, D. S., Kleinen, T., Lambert, F., van der Leeuw, S., Liddy, H., Loutre, M.F., McGee, D., Rehfeld, K., Rhodes, R., Seddon, A. W. R., Trauth, M. H., Vanderveken, L., and Yu, Z.: Past abrupt changes, tipping points and cascading impacts in the Earth system, Nat. Geosci., 14, 550–558, 2021. a, b
Brown, N. and Galbraith, E. D.: Hosed vs. unhosed: interruptions of the Atlantic Meridional Overturning Circulation in a global coupled model, with and without freshwater forcing, Clim. Past, 12, 1663–1679, https://doi.org/10.5194/cp1216632016, 2016. a
Bury, T. M., Bauch, C. T., and Anand, M.: Detecting and distinguishing tipping points using spectral early warning signals, J. Roy. Soc. Interface, 17, 20200482, https://doi.org/10.1098/rsif.2020.0482, 2020. a, b, c, d, e, f
Caesar, L., Rahmstorf, S., Robinson, A., Feulner, G., and Saba, V.: Observed fingerprint of a weakening Atlantic Ocean overturning circulation, Nature, 556, 191–196, 2018. a
Capron, E., Landais, A., Chappellaz, J., Schilt, A., Buiron, D., DahlJensen, D., Johnsen, S. J., Jouzel, J., LemieuxDudon, B., Loulergue, L., Leuenberger, M., MassonDelmotte, V., Meyer, H., Oerter, H., and Stenni, B.: Millennial and submillennial scale climatic variations recorded in polar ice cores over the last glacial period, Clim. Past, 6, 345–365, https://doi.org/10.5194/cp63452010, 2010. a, b, c, d, e, f
Carpenter, S. R. and Brock, W. A.: Rising variance: a leading indicator of ecological transition, Ecol. Lett., 9, 311–318, 2006. a
Centre for Ice and Climate: Data, icesamples and software, Københavns Universitet, Centre for Ice and Climate, Niels Bohr Institute [data set], https://www.iceandclimate.nbi.ku.dk/data/ (last access: 20 March 2024), 2024. a
Cimatoribus, A. A., Drijfhout, S. S., Livina, V., and van der Schrier, G.: Dansgaard–Oeschger events: bifurcation points in the climate system, Clim. Past, 9, 323–333, https://doi.org/10.5194/cp93232013, 2013. a
Clements, C. F. and Ozgul, A.: Rate of forcing and the forecastability of critical transitions, Ecol. Evol., 6, 7787–7793, 2016. a
Cleveland, W., Grosse, E., and Shyu, W.: Local regression models. Chapter 8 in Statistical models in S (JM Chambers and TJ Hastie eds.), 608 p., Wadsworth & Brooks/Cole, Pacific Grove, CA, 1992. a
Dakos, V., Scheffer, M., van Nes, E. H., Brovkin, V., Petoukhov, V., and Held, H.: Slowing down as an early warning signal for abrupt climate change, P. Natl. Acad. Sci. USA, 105, 14308–14312, 2008. a, b, c
Dakos, V., Carpenter, S. R., Brock, W. A., Ellison, A. M., Guttal, V., Ives, A. R., Kéfi, S., Livina, V., Seekell, D. A., van Nes, E. H., and Scheffer, M.: Methods for detecting early warnings of critical transitions in time series illustrated using simulated ecological data, PloS one, 7, e41010, https://doi.org/10.1371/journal.pone.0041010, 2012. a, b, c, d, e, f
Dansgaard, W., Johnsen, S., Clausen, H., DahlJensen, D., Gundestrup, N., Hammer, C., Hvidberg, C., Steffensen, J., Sveinbjörnsdottir, A., Jouzel, J., and Bond, G.: Evidence for general instability of past climate from a 250kyr icecore record, Nature, 364, 218–220, 1993. a, b
Desroches, M., Guckenheimer, J., Krauskopf, B., Kuehn, C., Osinga, H. M., and Wechselberger, M.: Mixedmode oscillations with multiple time scales, Siam Rev., 54, 211–288, 2012. a, b
Ditlevsen, P. and Ditlevsen, S.: Warning of a forthcoming collapse of the Atlantic meridional overturning circulation, Nat. Commun., 14, 1–12, 2023. a, b
Ditlevsen, P. D. and Johnsen, S. J.: Tipping points: early warning and wishful thinking, Geophys. Res. Lett., 37, L19703, https://doi.org/10.1029/2010GL044486, 2010. a, b
Ditlevsen, P. D., Svensmark, H., and Johnsen, S.: Contrasting atmospheric and climate dynamics of the lastglacial and Holocene periods, Nature, 379, 810, https://doi.org/10.1038/379810a0, 1996. a
Dokken, T. M., Nisancioglu, K. H., Li, C., Battisti, D. S., and Kissel, C.: DansgaardOeschger cycles: Interactions between ocean and sea ice intrinsic to the Nordic seas, Paleoceanography, 28, 491–502, 2013. a, b
EPICA community members: High resolution record of Northern Hemisphere climate extending into the last interglacial period, Nature, 431, 147–151, 2004. a, b
FitzHugh, R.: Impulses and physiological states in theoretical models of nerve membrane, Biophys. J., 1, 445, https://doi.org/10.1016/S00063495(61)869026, 1961. a, b
Gkinis, V., Simonsen, S. B., Buchardt, S. L., White, J., and Vinther, B. M.: Water isotope diffusion rates from the NorthGRIP ice core for the last 16,000 years–Glaciological and paleoclimatic implications, Earth Planet. Sc. Lett., 405, 132–141, 2014. a, b
Held, H. and Kleinen, T.: Detection of climate system bifurcations by degenerate fingerprinting, Geophys. Res. Lett., 31, L23207, https://doi.org/10.1029/2004GL020972, 2004. a
Henry, L., McManus, J. F., Curry, W. B., Roberts, N. L., Piotrowski, A. M., and Keigwin, L. D.: North Atlantic ocean circulation and abrupt climate change during the last glaciation, Science, 353, 470–474, 2016. a
IPCC: Framing, Context, and Methods. In Climate Change 2021 – The Physical Science Basis: Working Group I Contribution to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, chap. 1, 147–286, Cambridge University Press, https://doi.org/10.1017/9781009157896.003, 2023. a
Kilbourne, K. H., Wanamaker, A. D., MoffaSanchez, P., Reynolds, D. J., Amrhein, D. E., Butler, P. G., Gebbie, G., Goes, M., Jansen, M. F., Little, C. M., Mette, M., MorenoChamarro, E., Ortega, P., OttoBliesner, B. L., Rossby, T., Scourse, J., and Whitney, N. M.: Atlantic circulation change still uncertain, Nat. Geosci., 15, 165–167, 2022. a
Kindler, P., Guillevic, M., Baumgartner, M., Schwander, J., Landais, A., and Leuenberger, M.: Temperature reconstruction from 10 to 120 kyr b2k from the NGRIP ice core, Clim. Past, 10, 887–902, https://doi.org/10.5194/cp108872014, 2014. a
Klockmann, M., Mikolajewicz, U., Kleppin, H., and Marotzke, J.: Coupling of the subpolar gyre and the overturning circulation during abrupt glacial climate transitions, Geophys. Res. Lett., 47, e2020GL090361, https://doi.org/10.1029/2020GL090361, 2020. a
Klus, A., Prange, M., Varma, V., Tremblay, L. B., and Schulz, M.: Abrupt cold events in the North Atlantic Ocean in a transient Holocene simulation, Clim. Past, 14, 1165–1178, https://doi.org/10.5194/cp1411652018, 2018. a
Koper, M. T.: Bifurcations of mixedmode oscillations in a threevariable autonomous Van der PolDuffing model with a crossshaped phase diagram, Physica D, 80, 72–94, 1995. a
Kuehn, C.: A mathematical framework for critical transitions: normal forms, variance and applications, J. Nonlin. Sci., 23, 457–510, 2013. a
Kuniyoshi, Y., AbeOuchi, A., SherriffTadano, S., Chan, W.L., and Saito, F.: Effect of Climatic Precession on DansgaardOeschgerLike Oscillations, Geophys. Res. Lett., 49, e2021GL095695, https://doi.org/10.1029/2021GL095695, 2022. a
Kwasniok, F.: Analysis and modelling of glacial climate transitions using simple dynamical systems, Philos. T. Roy. Soc. Lond. A, 371, 20110472, https://doi.org/10.1098/rsta.2012.0374, 2013. a, b
Lenton, T. M., Livina, V. N., Dakos, V., and Scheffer, M.: Climate bifurcation during the last deglaciation?, Clim. Past, 8, 1127–1139, https://doi.org/10.5194/cp811272012, 2012. a
Li, C. and Born, A.: Coupled atmosphereiceocean dynamics in DansgaardOeschger events, Quaternary Sci. Rev., 203, 1–20, 2019. a
Livina, V. N. and Lenton, T. M.: A modified method for detecting incipient bifurcations in a dynamical system, Geophys. Res. Lett., 34, L03712, https://doi.org/10.1029/2006GL028672, 2007. a
Lohmann, J. and Ditlevsen, P. D.: A consistent statistical model selection for abrupt glacial climate changes, Clim. Dynam., 52, 6411–6426, 2019. a, b
Lohmann, J. and Ditlevsen, P. D.: Risk of tipping the overturning circulation due to increasing rates of ice melt, P. Natl. Acad. Sci. USA, 118, e2017989118, https://doi.org/10.1073/pnas.2017989118, 2021. a
Lohmann, J., Dijkstra, H. A., Jochum, M., Lucarini, V., and Ditlevsen, P. D.: Multistability and Intermediate Tipping of the Atlantic Ocean Circulation, arXiv preprint arXiv:2304.05664, 2023. a
Lucarini, V. and Stone, P. H.: Thermohaline circulation stability: A box model study. Part I: Uncoupled model, J. Climate, 18, 501–513, 2005. a
MalmiercaVallet, I., Sime, L. C., and the D–O community members: Dansgaard–Oeschger events in climate models: review and baseline Marine Isotope Stage 3 (MIS3) protocol, Clim. Past, 19, 915–942, https://doi.org/10.5194/cp199152023, 2023. a
Martrat, B., Grimalt, J. O., LopezMartinez, C., Cacho, I., Sierro, F. J., Flores, J. A., Zahn, R., Canals, M., Curtis, J. H., and Hodell, D. A.: Abrupt temperature changes in the Western Mediterranean over the past 250,000 years, Science, 306, 1762–1765, 2004. a
MassonDelmotte, V., Zhai, P., Pirani, A., Connors, S. L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M. I., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J. B. R., Maycock, T. K., Waterfield, T., Yelekci, O., Yu, R., and Zhou, B.: Climate change 2021: the physical science basis, Contribution of working group I to the sixth assessment report of the intergovernmental panel on climate change, Cambridge University Press, Cambridge, UK and New York, NY, USA, https://doi.org/10.1017/9781009157896, in press, 2021. a
Meisel, C. and Kuehn, C.: Scaling effects and spatiotemporal multilevel dynamics in epileptic seizures, PLoS One, 7, e30371, https://doi.org/10.1371/journal.pone.0030371, 2012. a, b
Menviel, L. C., Skinner, L. C., Tarasov, L., and Tzedakis, P. C.: An ice–climate oscillatory framework for Dansgaard–Oeschger cycles, Nat. Rev. Earth Environ., 1, 677–693, 2020. a
Michel, S. L., Swingedouw, D., Ortega, P., Gastineau, G., Mignot, J., McCarthy, G., and Khodri, M.: Early warning signal for a tipping point suggested by a millennial Atlantic Multidecadal Variability reconstruction, Nat. Commun., 13, 5176, https://doi.org/10.1038/s41467022327043, 2022. a, b
Mitsui, T.: takahito321/PredictabilityofDOcooling: Release, Zenodo [code], https://doi.org/10.5281/zenodo.10841655, 2024. a
Mitsui, T. and Crucifix, M.: Influence of external forcings on abrupt millennialscale climate changes: a statistical modelling study, Clim. Dynam., 48, 2729–2749, 2017. a, b
Nagumo, J., Arimoto, S., and Yoshizawa, S.: An active pulse transmission line simulating nerve axon, Proc. IRE, 50, 2061–2070, 1962. a, b
O'Sullivan, E., Mulchrone, K., and Wieczorek, S.: Rateinduced tipping to metastable zombie fires, P. Roy. Soc. A, 479, 20220647, https://doi.org/10.1098/rspa.2022.0647, 2023. a
Peltier, W. R. and Vettoretti, G.: DansgaardOeschger oscillations predicted in a comprehensive model of glacial climate: A “kicked” salt oscillator in the Atlantic, Geophys. Res. Lett., 41, 7306–7313, 2014. a, b
Rahmstorf, S.: Ocean circulation and climate during the past 120,000 years, Nature, 419, 207–214, 2002. a, b
Rahmstorf, S., Box, J. E., Feulner, G., Mann, M. E., Robinson, A., Rutherford, S., and Schaffernicht, E. J.: Exceptional twentiethcentury slowdown in Atlantic Ocean overturning circulation, Nat. Clim. Change, 5, 475–480, 2015. a
Rasmussen, S. O., Abbott, P. M., Blunier, T., Bourne, A. J., Brook, E., Buchardt, S. L., Buizert, C., Chappellaz, J., Clausen, H. B., Cook, E., DahlJensen, D., Davies, S. M., Guillevic, M., Kipfstuhl, S., Laepple, T., Seierstad, I. K., Severinghaus, J. P., Steffensen, J. P., Stowasser, C., Svensson, A., Vallelonga, P., Vinther, B. M., Wilhelms, F., and Winstrup, M.: A stratigraphic framework for abrupt climatic changes during the Last Glacial period based on three synchronized Greenland icecore records: refining and extending the INTIMATE event stratigraphy, Quaternary Sci. Rev., 106, 14–28, 2014. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o
Rial, J. and Yang, M.: Is the frequency of abrupt climate change modulated by the orbital insolation?, Geophysical monographAmerican Geophysical Union, 173, 167–174, 2007. a, b
Riechers, K., Mitsui, T., Boers, N., and Ghil, M.: Orbital insolation variations, intrinsic climate variability, and Quaternary glaciations, Clim. Past, 18, 863–893, https://doi.org/10.5194/cp188632022, 2022. a, b
Ritchie, P. D. L., Alkhayuon, H., Cox, P. M., and Wieczorek, S.: Rateinduced tipping in natural and human systems, Earth Syst. Dynam., 14, 669–683, https://doi.org/10.5194/esd146692023, 2023. a
Roberts, A. and Saha, R.: Relaxation oscillations in an idealized ocean circulation model, Clim. Dynam., 48, 2123–2134, 2017. a, b
Ruth, U., Wagenbach, D., Steffensen, J. P., and Bigler, M.: Continuous record of microparticle concentration and size distribution in the central Greenland NGRIP ice core during the last glacial period, J. Geophys. Res.Atmos., 108, 4091, https://doi.org/10.1029/2002JD002376, 2003. a, b
Rypdal, M.: Earlywarning signals for the onsets of Greenland interstadials and the Younger Dryas–Preboreal transition, J. Climate, 29, 4047–4056, 2016. a, b, c
Sadatzki, H., Dokken, T. M., Berben, S. M., Muschitiello, F., Stein, R., Fahl, K., Menviel, L., Timmermann, A., and Jansen, E.: Sea ice variability in the southern Norwegian Sea during glacial DansgaardOeschger climate cycles, Science Adv., 5, eaau6174, https://doi.org/10.1126/sciadv.aau6174, 2019. a
Sakai, K. and Peltier, W. R.: A dynamical systems model of the DansgaardOeschger oscillation and the origin of the Bond cycle, J. Climate, 12, 2238–2255, 1999. a
Scheffer, M., Bascompte, J., Brock, W. A., Brovkin, V., Carpenter, S. R., Dakos, V., Held, H., Van Nes, E. H., Rietkerk, M., and Sugihara, G.: Earlywarning signals for critical transitions, Nature, 461, 53–59, 2009. a, b, c, d
Seierstad, I. K., Abbott, P. M., Bigler, M., Blunier, T., Bourne, A. J., Brook, E., Buchardt, S. L., Buizert, C., Clausen, H. B., Cook, E., DahlJensen, D., Davies, S., Guillevic, M., Johnsen, S. J., Pedersen, D. S., Popp, T. J., Rasmussen, S. O., Severinghaus, J., Svensson, A., and Vinther, B. M.: Consistently dated records from the Greenland GRIP, GISP2 and NGRIP ice cores for the past 104 ka reveal regional millennialscale δ^{18}O gradients with possible Heinrich event imprint, Quaternary Sci. Rev., 106, 29–46, 2014. a, b, c, d, e
Stommel, H.: Thermohaline convection with two stable regimes of flow, Tellus, 13, 224–230, 1961. a, b
Strogatz, S. H. (Ed.): Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering, CRC Press, https://doi.org/10.1201/9780429399640, 2018. a, b
Theiler, J., Eubank, S., Longtin, A., Galdrikian, B., and Farmer, J. D.: Testing for nonlinearity in time series: the method of surrogate data, Physica D, 58, 77–94, 1992. a, b
Thomas, Z. A., Kwasniok, F., Boulton, C. A., Cox, P. M., Jones, R. T., Lenton, T. M., and Turney, C. S. M.: Early warnings and missed alarms for abrupt monsoon transitions, Clim. Past, 11, 1621–1633, https://doi.org/10.5194/cp1116212015, 2015. a
Thompson, J. M. T. and Sieber, J.: Climate tipping as a noisy bifurcation: a predictive technique, IMA J. Appl. Math., 76, 27–46, 2011. a, b
Timmermann, A., Gildor, H., Schulz, M., and Tziperman, E.: Coherent resonant millennialscale climate oscillations triggered by massive meltwater pulses, J. Climate, 16, 2569–2585, 2003. a
van der Bolt, B., van Nes, E. H., and Scheffer, M.: No warning for slow transitions, J. Roy. Soc. Interface, 18, 20200935, https://doi.org/10.1098/rsif.2020.0935, 2021. a
Vettoretti, G. and Peltier, W. R.: Fast physics and slow physics in the nonlinear Dansgaard–Oeschger relaxation oscillation, J. Climate, 31, 3423–3449, 2018. a, b
Vettoretti, G., Ditlevsen, P., Jochum, M., and Rasmussen, S. O.: Atmospheric CO2 control of spontaneous millennialscale ice age climate oscillations, Nat. Geosci., 15, 300–306, 2022. a, b, c
Wieczorek, S., Xie, C., and Ashwin, P.: Rateinduced tipping: Thresholds, edge states and connecting orbits, Nonlinearity, 36, 3238, https://doi.org/10.1088/13616544/accb37, 2023. a
Yiou, R., Fuher, K., Meeker, L., Jouzel, J., Johnsen, S., and Mayewski, P. A.: Paleoclimatic variability inferred from the spectral analysis of Greenland and Antarctic icecore data, J. Geophys. Res.Oceans, 102, 26–441, 1997. a
Zhang, X., Barker, S., Knorr, G., Lohmann, G., Drysdale, R., Sun, Y., Hodell, D., and Chen, F.: Direct astronomical influence on abrupt climate variability, Nat. Geosci., 14, 819–826, 2021. a
 Abstract
 Introduction
 Data and methods
 Results
 Possible dynamical mechanisms
 Summary and discussion
 Appendix A: Probability of observing robust precursor signals
 Appendix B: Details of conceptual models used in Fig. 5
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Data and methods
 Results
 Possible dynamical mechanisms
 Summary and discussion
 Appendix A: Probability of observing robust precursor signals
 Appendix B: Details of conceptual models used in Fig. 5
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement