Dansgaard-Oeschger events: tipping points in the climate system

Dansgaard-Oeschger events are a prominent mode of variability in the records of the last glacial cycle. Various prototype models have been proposed to explain these rapid climate fluctuations, and no agreement has emerged on which may be the more correct for describing the paleoclimatic signal. In this work, we assess the bimodality of the system reconstructing the topology of the multi--dimensional attractor over which the climate system evolves. We use high-resolution ice core isotope data to investigate the statistical properties of the climate fluctuations in the period before the onset of the abrupt change. We show that Dansgaard-Oeschger events have weak early warning signals if the ensemble of events is considered. We find that the statistics are consistent with the switches between two different climate equilibrium states in response to a changing external forcing (e.g. solar, ice sheets...), either forcing directly the transition or pacing it through stochastic resonance. These findings are most consistent with a model that associates Dansgaard-Oeschger with changing boundary conditions, and with the presence of a bifurcation point.


I. INTRODUCTION
The largest variability in temperature over the last sixty thousand years is connected with Dansgaard-Oeschger events (DOs) [1,2]. These are fast warming episodes (in the North Atlantic region 5 ÷ 10 • C in a few decades), followed by a gradual cooling that lasts from hundreds to thousands of years, often with a final jump back to stadial condition. They occurred with a periodicity of approximately 1500 yr [3]. The relation between DOs and large changes in the Atlantic meridional overturning circulation is generally considered well established, despite the limitations of the paleoclimatic records [3][4][5][6][7]. Various simple models have been proposed to explain these rapid climate fluctuations [3,6,7], but until now no observational constraint has been forwarded to choose between different theories. Here, we use high-resolution ice core isotope data [8,9] to investigate the statistical properties of the climate fluctuations [10][11][12] in the period before the onset of the abrupt change. We analyse δ 18 O isotope data from the NGRIP ice core [8,9]. The data spans the time interval from 59 420 to 14 777 yr before 2000 A.D. (years b2k) with an average resolution of 2.7 yr (more than 80 % of the data has a resolution better than 3.5 yr). δ 18 O can be interpreted as a proxy for atmospheric temperature over Greenland. We use the dating for the onset of the interstadials given in Svensson et al. [9]. The dataset then spans the events numbered 2 to 16. In the first part of the paper, we discuss the bimodality of the time series. In particular, we establish that bimodality is a robust feature of the climate system that produced it, and not an artifact of the projection of complex dynamics onto a scalar quantity. This is done using a phase embedding technique well known in the non-linear dynamics community, but not often used in climate sciences.
In the second part of the paper, we show that the statistical properties of the paleoclimatic time series can be used to distinguish between the various models that have been proposed to explain DOs. Despite limitations of the available data, we find that the statistics are most compatible with a system that switches between two different climate equilibrium states in response to a changing external forcing. Other hypotheses [noise-induced transitions and autonomous oscillations, suggested by 3, 6, 7, 13] appear less compatible with the data. The study is also meant to demonstrate that the average statistical properties of fluctuations of a time series can be used to extract precious information on the dynamics of the system that produced it.

II. BIMODALITY
Often, bimodality is implicitly assumed for Dansgaard-Oeschger events, i.e. it is assumed that the climate system can switch between two distinct equilibrium states. The bimodality hypothesis for this climate proxy has been tested in previous works [14][15][16], but these make the strong assumption that the bimodality of the recorded time series can be directly ascribed to the climate system that produced it. This is not necessarily the case if the evolution of a system with many degrees of freedom, as the real climate system, is projected on a scalar record. Exploiting Takens' embedding theorem [17] for reconstructing the phase space (the space that includes all the states of a system), we assess bimodality without making any assumption on the underlying model beforehand. This approach has the advantage of avoiding artifacts that may arise from the projection of a higher dimensional system onto a scalar time series.

A. Phase space reconstruction technique
The technique that we apply is outlined in Abarbanel [18]. First, the irregular time series is linearly interpolated to a regular time step, 4 yr in the results shown here. Data after the last DO is discarded for producing the phase space reconstruction. After detrending, an optimal time lag has to be empirically chosen, guided by the first minimum of the average mutual information function (20 yr in our case). With copies of the regular time series lagged by multiples of the optimal time lag as coordinates, a phase space reconstruction is obtained. The dimension of the embedding space is, in general, greater than or equal to the one of the physical system which produced the signal, and the dynamics of the system are thus completely unfolded. The number of dimensions is chosen as the minimum one that brings the fraction of empirical global false neighbours to the background value [19]. In our reconstructions, four dimensions are used. We define as global false neighbours in dimension n those couples of points whose distance increases more than 15 times, or more than two times the standard deviation of the data, when the embedding dimension is increased by one. The Scientific Tools for Python (SciPy) implementation of the kd-tree algorithm is used for computing distances in the n-dimensional space. The 4-dimensional phase space reconstruction is presented in Fig. 1, performed separately for the data before and after year 22 000 b2k. The results are converted to a polar coordinate system and the anomaly from the average is considered, following Stephenson et al. [20]. Subsequently, the PDFs for the coordinates distributions are computed through a gaussian kernel estimation (blue line in Fig. 1). The error of the PDFs is defined as the standard deviation of an ensemble of 50 000 synthetic datasets (shaded region in the figure), obtained by bootstrapping [see e.g. 21]. The theoretical distributions (black lines) refer to a multidimensional normal distribution (the simplest null hypothesis, 20).

B. Results: bimodality and bifurcation
The δ 18 O isotope data possess a distinct trend that marks the drift of the glacial climate towards the last glacial maximum (approximately 22 000 yr ago), see Fig. 2 (top). Also, associated with the drift is a gradual increase in variance that makes up an early warning signal (EWS) for the abrupt change towards the last glacial termination [approximately 18 000 yr ago, 12,22]. It turns out that the trend enhances the bimodality of the data but detrending is required to unequivocally attribute this bimodality to the DOs, as the trend is connected with a longer time scale process.
In the four left panels of Fig. 1 (before year 22 000 b2k), the normal distribution clearly lies outside the error bars for all four dimensions, and has to be rejected. The bimodality of the data is especially evident from the PDFs of the angular coordinates. The radial distribution has a heavy tail which suggests a concentration of data far from the main peak. These features are very robust, and they can be unambiguously detected also in phase space reconstructions using different interpolation steps, time lags and number of dimensions. In the four right panels (after year 22 000 b2k), the deviations from normality are much smaller, and become indistinguishable from a normal distribution when the errors are taken into account. The only statistically non insignificant deviations from normality after year 22 000 b2k are seen in panel h, which represents the longitude coordinate of the polar system. However, no bimodality is observed, but only a non spherically-symmetrical distribution of the data. On this basis, we can confirm the claim of Livina et al. [14,15]; around year 22 000 b2k undergoes a transition from a bimodal to a unimodal behaviour, possibly linked to a bifurcation.

III. EXTERNAL FORCING, NOISE OR AUTONOMOUS DYNAMICS?
We want now to analyse the statistical properties of the shorter time-scale fluctuations in the ice core record that occur before the onset of DOs. This aim is more ambitious than the assessment of bimodality, and pushes to the extreme the amount of information that can be extracted from the data. Despite the limits determined by the quality of the paleoclimatic time series, the analysis allows to falsify various prototype models that have been proposed as an explanation for DOs, while at the same time it suggests a mechanism that is compatible with the statistical properties of the time series. Such a selection is possible since the different prototype models have different, peculiar, "noise signatures" before an abrupt transition takes place. It must be stressed that the problem is an inverse one; multiple mathematical models qualitatively consistent with the "noise signatures" can be developed, and only physical mechanisms can guide the final choice.
In several recent works [10][11][12][13], an abrupt transition is treated as the jump between different steady states in a system that can be decomposed in a low dimensional non-linear component and a small additive white noise component. The jump is considered to be due to the approach to a bifurcation point, usually a fold bifurcation [see e.g. 23, for a discussion of fold bifurcations], in which a steady state loses its stability in response to a change in the boundary conditions. If a system perturbed by a small amount of additive white noise approaches a fold bifurcation point, an increase in the variance is observed, as well as in its autocorrelation [10,12,13,24]. In addition, Detrended Fluctuation Analysis (DFA) can be used to detect the approach to a bifurcation point [11]. This tool, equivalent to others that measure the decay rate of fluctuations in a nonlinear system [25,26], has the clear advantage over c and σ 2 of inherently distinguishing between the actual fluctuations of the system and trends that may be part of the signal. This follows from the fact that the detrending is part of the DFA procedure itself, and different orders of detrending can be tested, provided that a sufficient amount of data is available to sustain the higher order fit. Approaching a bifurcation point, c and α tend to increase towards 1 and 3/2 respectively [10,11]. This reflects the fact that the fluctuations in the system are increasingly correlated as the bifurcation point is approached: they become more similar to a Brownian motion as the restoring force of the system approaches zero.
We thus focus on three characteristics of the data: the correlation coefficient of the time series at lag-1 (c), the variance (σ 2 ), and the DFA exponent (α) [26]. These quantities are related to the decay rate of the shorter time-scale fluctuations, and their average magnitude. A similar type of analysis has been used before in the search of EWSs for abrupt climate transitions [see e.g. [10][11][12][13]. Here we use EWSs for the first time to put observational constraints to the various prototype models for DOs. In this way, we can select which prototype model is compatible with the EWSs that are contained in the data. In contrast to earlier studies we consider an EWS not as a precursor of a single event, but analyse instead the average characteristics of the EWS found in the ensemble of DOs. This approach is much more robust, since EWSs have a predictable behaviour only when computed as an average over an ensemble of events, as demonstrated by Kuehn [24].

A. Prototype models and early warning signals
The reason why we can discriminate between various prototype models is that, in the case of a system that crosses, or comes close to a tipping point, σ 2 , c and α all increase before abrupt change occurs, while for other prototype models this is not the case.
For each prototype model that we considered, the relevant equation, or set of equations, is integrated with an Euler-Maruyama method [27]. The associated time series is then analysed with respect to lag-1 autocorrelation, variance, and the DFA exponent (details of the analyses are given in Sect. III B).

Tipping point in a double well
EWSs are visible if the DOs are caused by a slowly changing external forcing, in which a critical parameter, like northward freshwater transport, varies over a large enough range for the meridional overturning circulation of the ocean to cross the two tipping points that mark the regime of multiple equilibria. In physical terms, this model implies that a changing external forcing is causing a previously stable state to disappear. In Fig. 3a this situation is sketched.
As a prototype for a system crossing a tipping point, the same equations used in [13] are used. In this case, the evolution of the system can be thought of as the motion of a sphere on a potential surface that has two possible equilibrium states, divided by a potential ridge. This potential surface can change in time, due to slow changes in the forcing. Beyond a critical point the system can jump to a different equilibrium (Fig. 4).

Noise induced transition
If the external forcing is constant in time, abrupt transitions can still occur if more than one equilibrium state is available, and noise is sufficiently strong to trigger the transition. An example of this case, in the simple case of a fixed double well with constant white noise, is shown in Fig. 5. In this model, the system is normally in the cold stadial state but can jump due to the noise in the forcing, say atmospheric weather, to a warmer, interstadial state. The change in external forcing is considered too small to bring the system close to a tipping point. In this case no trend is visible in the statistical properties of the time series, thus no EWS can be detected. This is the view of Ganopolski and Rahmstorf [3], Ditlevsen and Johnsen [13].

Other models
We will show that other models proposed in the literature for explaining DOs do not posses any EWS. No EWS can be detected if DOs are due to an unforced (ocean) mode of oscillation [7,28] since, for an autonomous oscillation, there are in general no changes in the stability properties of the climate system while it evolves in time [see e.g. 29]. DOs have been linked to oscillations close to a "homoclinic orbit" by Timmermann et al. [6]. A homoclinic orbit is a periodic solution to a system of Ordinary Differential Equations that has infinite period, but in the presence of a small amount of noise it features finite return times. A sketch of this prototype model is shown in Fig. 3b. To describe this model the set of Ordinary Differential Equations given in Crommelin et al. [30] is used. This minimal mathematical model has been suggested before in the context of atmosphere regime behaviour Charney and DeVore [31], but it is mathematically equivalent to the mechanism for DOs suggested in Timmermann et al. [6], Abshagen and Timmermann [32].
The system of ordinary differential equations (2.6) in [30] is integrated to investigate EWSs for a periodic cycle close to a homoclinic orbit (Fig. 6). This model is chosen as it provides a simple description of the results of [6], having a very similar bifurcation diagram (the fundamental element is the presence of a fold-Hopf bifurcation). The parameter choice is the same as in [30], with x * 1 = 0.89 and r = −0.925. As the time series is strongly non-stationary, linear DFA is not a good measure of the fluctuation decay in this case, and higher order detrending is needed (orders higher than 2 give results consistent with the case of quadratic detrending).
After integrating these equations, again no precursor for abrupt changes can be detected (Fig. 6). It must be stressed that, in order to correctly compute the statistical properties of the time series, it is of fundamental importance that the data are detrended (linear or higher order detrending), to ensure that the findings can be ascribed to random fluctuations and not to the non-stationarity of the data. In particular, in the case of the DFA exponent computation for the latter prototype model, the strong non-stationarity of the data requires a quadratic detrending.
The results from other models that have been suggested as a prototype for oscillations close to a homoclinic orbit [33,34] consistently show no EWS before abrupt transitions (Fig. 7). The Welander oscillator [33] was suggested as a prototype model for DOs in Abshagen and Timmermann [32]. We considered Eq. (4) of Abshagen and Timmermann [32] with µ = 0.73, ν 1 = 0.1, ν 2 = 3.0 and σ = 8. In the equation, ν is defined as a Heavyside function as this choice gives a time series qualitatively closer to the paleoclimatic record. Also in this case, no EWSs are detected (Fig. 7).

B. Fluctuations analysis
We now have demonstrated that EWSs are characteristic features of only one prototype model. If the ice core record also features EWSs, the other prototype models are falsified and we can attribute the DOs to the mechanism which implies the crossing of a tipping point. The ice core data are shown in Fig. 2. The EWSs of the time series are computed in a sliding window of 250 yr wide, and the results are plotted at the right end of the sliding window to avoid contamination from points after the DO onset. This window size corresponds on average to 100 data points for the variance computation. For c and α, the time series is first linearly interpolated to a step of 4 yr and approximately 60 points are used in the window. The time step is chosen in order to avoid overfitting the data even in the parts of the time series with lower resolution. The window width choice is a compromise between the need to minimise noise (larger window) and to maximise the number of independent points for the computation of trends (smaller window). The data is linearly detrended in the window before the computation, to compensate for potential nonstationarities in the time series. In the computation of the DFA exponent [26], 10 boxes lengths were used, ranging from 16 to 52 yr. The number of time lengths is limited for computational reasons. Different time steps for the interpolation and boxes lengths were tested. The results show small sensitivity to these choices, as long as the time step is kept below approximately 10 yr.
If one considers the successive DOs one by one (Fig. 2), it is clear that for some of the DOs a marked increase is found (e.g. the one at approximately 45 000 yr b2k), implying an EWS, but for other DOs no EWSs can be detected (e.g. the last one). This has to be expected from a stochastic system for which only the average behaviour is predictable Kuehn [24]. In other words, there is a strong element of randomness associated with DOs; the external forcing that changes the critical parameters of the system may be non-periodic, and also the noise associated with shorter term weather and climate fluctuations introduces an element of unpredictability.

C. Results
As discussed by Kuehn [24], analysing the properties of an ensemble of events, instead of a single realisation, is an approach more robust and better justified on theoretical grounds. We thus use EWSs to characterise the property of the "mean DO" instead of trying to predict the onset of the following transition. For this reason we consider the whole ensemble of DOs, instead of each single DO. To this end we transform the time axis to perform a more robust analysis on the data. The time series is now cut into slices that end 100 yr after the transition onset and start either 100 yr after the previous DO onset or, if the time span between two DOs is sufficiently long, 2900 yr before the following DO onset (Fig. 8a). The onset of the transitions then are translated to year −100. In this way, an ensemble of DOs is defined, being composed of slices of the original time series of variable length, synchronised at the onset of each DO. If the quantities used for EWS detection are then averaged for the whole ensemble, a moderate but significant increase is observed in all three fluctuation properties, starting approximately from year −1800 until year −250 (Fig. 8b-d). The standard deviation of the ensemble is large, in particular at times far from the DO onset because the ensemble has fewer members there, but the trend can not be discarded as a random fluctuation for all three cases. To prove this, a linear least square fit of the data in the interval −1800 to −250 is performed, only using independent points from the ensemble (farther than 250 yr apart), and obtaining an error interval from a bootstrapped ensemble of 50 000 members. The results of this fitting procedure are reported in Table I. In all three cases the linear trends are well above their standard deviation, providing a strong indication of the robustness of the signal. These results are consistent with a scenario in which DOs are associated with the crossing of a tipping point, in contrast with suggestions from previous work [13]. In order to check the robustness of the findings on the order of detrending for DFA, the computation has been repeated with a quadratic detrending, actually giving a slightly stronger trend (see Table I).
The fact that c and α do not reach their theoretical values for a critical point, respectively 1 and 3/2, can easily be explained: the noise in the system triggers the transition before the actual tipping point is reached [35]. Also, the three quantities c, σ 2 and α decrease again in the last few hundred years before the DO onset. To explain this decrease, we refer to Kuehn [24], who shows that this behaviour can be expected from a system approaching a bifurcation point. For several idealised systems containing different forms of tipping points (included the fold bifurcation considered here), a decrease in variance in the immediate vicinity of the bifurcation point is found, while the well known increase is observed when the bifurcation point is farther away in the parameter space. The observation of a decrease in variance just before the DO onset, after a longer period of increase, is an important lesson for the practical application of these techniques in the detection of approaching tipping points.

IV. SUMMARY
In this work, we performed two sets of analysis on a well-known paleoclimatic record, the δ 18 O isotope data from the NGRIP ice core [8,9], believed to be representative of temperature over the North Atlantic sector. We assessed bimodality of the system using a phase-space embedding technique, that guarantees that the bimodality is not an artifact of projection of the complex dynamics of the climate system on a scalar time series. We confirm with this technique the claim of Livina et al. [14,15], that a switch from bimodal to unimodal behaviour is observed in the time series around year 22 000 before present.
Secondly, we analysed the statistical properties of the fluctuations in the paleoclimatic time series, before the onset of DOs. In particular, we focused on the average properties of the events considered as an ensemble instead of each separately. Despite the high level of noise in the data, we were able to conclude that EWSs are present, consistently with the hypothesis that DOs result from a repeated crossing of a tipping point in the climate system. Other prototype models that have been proposed [3,6,7,13] are less consistent with the data, as their mechanisms do not involve any transition from which EWSs can, at least in general, be expected. Ditlevsen and Johnsen [13] came to opposite conclusions, but they did not considered the average behaviour of the ensemble, while we think that this may be a step of fundamental importance.
A disclaimer has to be made: our conclusions hold for the ensemble of the events, and are a probabilistic statement: we can only claim that a scenario that does not include a tipping point with EWS is unlikely. The trends of Table I are between two and three standard deviations apart from zero. This means that the probability that the real trend in each indicator is zero is low but non-zero, in the order of 1 %, assuming a normal distribution. Furthermore, given the rather complex techniques used, we can not rule out the possibility that the error estimates given in the Table  may be underestimated. A connection with thermohaline circulation instability remains in the domain of speculation, but is a plausible hypothesis considering the large evidence linking DOs signals to rapid changes in the meridional overturning circulation. Further investigation is needed to confirm this hypothesis and, more importantly, to address the fundamental question that remains open: if DOs are due to an external forcing, what is the nature of this forcing?
Looking beyond the results of the particular time-series used here, we suggest that EWSs may provide a useful guide for discriminating between different models meant to describe a climate signal. When data is scarce, the analysis of the average properties of the fluctuations can give important hints on the nature of the signal which produced it.  b-c and f-g are the PDFs for the two "latitude-like" coordinates. d and e are the PDF for the "longitude-like" coordinate. The black lines are the theoretical PDFs for a unimodal multinormal sample: in the polar reference system, a χ 2 distribution, with number of degrees of freedom equal to the number of dimensions, for the square of the radius (ρ 2 ). For the angular coordinates, a unimodal multinormal PDF would be proportional to sin 2 (φ1) and sin(φ2) respectively, while it would be uniform in φ3.  a Tipping points hypothesis. The system crosses the two tipping points (green triangles) periodically, due to an external forcing. The DOs give EWSs before the abrupt transition occurs. Overshooting is possible when overturning circulation recovers. b Homoclinic orbit model Timmermann et al. [6]. Here, DOs are due to the motion of the climate system close to a homoclinic orbit (grey) that connects to a periodic oscillation (Hopf bifurcation) Timmermann et al. [6], Crommelin et al. [30], Abshagen and Timmermann [32] (green circle). Tipping points are present in the climate system, but they do not determine the abrupt transitions. The climate is most of the time in stadial conditions. A possible path followed during the oscillation is shown in black for each case.  Fig. 4, but for the Welander oscillator. In this case, similarly to the one of the Charney and Devore model, the abrupt changes are an autonomous mode of the system. The first dimension of the 2-dimensional system is used. Also in this case a detrending of order higher than one improves the results from DFA, that ore otherwise strongly affected by the non-stationarity of the time series.  [9]. The time series are aligned so that the DOs start at year −100 (grey shading after event onset). Colours mark different events, from number 2 (blue) to number 16 (red). The time window used in the computations is drawn in black. Below, the lag-1 autocorrelation (b), variance (c) and linear detrended fluctuation analysis exponent (d). The blue line is the ensemble average, the shaded area marks the standard deviation. The dashed red line shows a least square regression of the data over the marked time range.