By comparing the high-resolution isotopic records from the Greenland Ice Core Project (GRIP) and the North Greenland Ice Core Project (NGRIP) ice cores, the common climate signal in the records has been approximately separated from local noise. From this, an objective criterion for defining Dansgaard–Oeschger (DO) events is achieved. The analysis identifies several additional short-lasting events, increasing the total number of DO events to 27 in the period 12–90 kyr before present (BP). The quasi-regular occurrence of the DO events could indicate a stochastic or coherent resonance mechanism governing their origin. From the distribution of waiting times, a statistical upper bound on the strength of a possible periodic forcing is obtained. This finding indicates that the climate shifts are purely noise driven with no underlying periodicity.
Abrupt temperature shifts between cold (stadial) states and warm (interstadial) states [Dansgaard–Oeschger (DO) events] are observed in the Greenland ice core isotope records (Dansgaard et al. 1993; Grootes et al. 1993). The transitions into the interstadials were abrupt and with temperature changes estimated from the paleo record of the order of 10°–15°C occurring within decades. These two distinct quasi-stable climate states are most likely linked to different modes of the Atlantic thermohaline circulation (Broecker et al. 1985; Stocker et al. 1992). The circulation in the warm state was similar to the present time circulation, while in the cold state the sinking took place at lower latitudes and to shallower depths (Alley and Clark 1999; Ganopolski and Rahmstorf 2001). A complete cessation of deep water formation is indicated by some models (Schmittner et al. 2002).
Identification of the mechanism causing these climate shifts is of primary importance for understanding the stability and mode of operation of this component of the climate system. Comparing the high-resolution isotope records from the Greenland Ice Core Project (GRIP) and the North Greenland Ice Core Project (NGRIP; North GRIP Members 2004, hereafter NGM), we can separate the climate signal from the local variability and the glaciological noise in the ice cores. We observe that for time scales longer than approximately 30 yr, the two ice cores are strongly correlated. Using the 30-yr-averaged record, we thus obtain an objective measure to define the DO events. The initiations and terminations of the DO events are defined from consecutive up-crossings and down-crossings through two levels in the anomaly record. The definition of the DO events by this procedure is quite robust with respect to the choice of anomaly levels. Comparing to the original visual numbering of the DO events (Dansgaard et al. 1993), we find several additional isotopic fluctuations that qualify as DO events. For example, DO event 2 consists of two closely spaced DO events, which based on the stratigraphy devised by Walker et al. (1999) and Björk et al. (1998) becomes Greenland Interstadial (GI) “GI2a” (youngest) and “GI2c” (oldest) separated by the stadial state “GI2b.” Using these conventions, the succeeding period is the Greenland Stadial (GS2). (For further explanation see the caption of Fig. 2.)
The isotope record from the Greenland Ice Sheet Project 2 (GISP2) ice core, dated using stratigraphic methods (Meese et al. 1997), shows a significant peak in the spectrum at 1470 yr (Yiou et al. 1997), indicating a possible periodic forcing of the climate system. Alley et al. (2001) proposed this to be a stochastic resonance, while Timmermann et al. (2003) proposed a coherent resonance phenomenon. The climate system itself is dominated by strongly fluctuating, irregular, fast time-scale noise, so it is highly unlikely that a strictly periodic signal would be internally generated. Only simplified models with few degrees of freedom exhibit strict cyclic behavior (Paillard and Labeyrie 1994). In a recent study, Roe and Steig (2004) noted that, 1470-yr periodicity aside, the Antarctic Byrd record is comparable with a simple autoregressive process, and the same is true for the Greenland records if an additional simple threshold rule is imposed.
A periodic component in the signal would likely be the result of a nonlinear response to a weak external periodic forcing, though the origin of such a forcing has not yet been identified. By observing the waiting times between consecutive DO events, it was noted that the record could be interpreted as having a preferred waiting time of 1470 yr and multiples of this period, corresponding to the system skipping a few transitions but still switching in phase with the external periodicity (Alley et al. 2001; Schulz 2002a; Rahmstorf 2003). Observing the waiting times rather than the power spectrum is advantageous in a noisy signal where the cyclicity is far from being sinusoidal. In this case, a large portion of the power will be in overtones, which possibly brings the power in the peak below detection into the noise level. This is the case for the ice core signal where the transitions into the interstadials are rapid and the transitions into the stadials are more gradual.
By comparing the ice core record with stochastic resonance models, we are able to estimate an upper bound on the strength of an external periodic forcing potential in comparison to the internal barrier between the two climate states. We find that with maximum likelihood the data are not a result of a process with a periodic component. Since there still remains discrepancy between the GRIP and the GISP2 datings, we have performed our analysis on both time scales. The GISP2 time scale agrees well with the uranium/thorium (238U/230Th) dating of the Chinese Hulu Cave stalagmite record (Wang et al. 2001), while the dating of the French Villars Cave stalagmite record (Genty et al. 2003) is somewhere in between the two ice core time scales. While we cannot reject the case of no periodic component, we can determine an upper limit of a periodic forcing such that a strength above this value can be rejected with 98% confidence. This upper limit varies slightly between the GRIP and the GISP2 time scales.
2. Defining the Dansgaard–Oeschger events
The power spectra of the GRIP and GISP2 records (Fig. 1) are estimated from the irregularly spaced data using the “Redfit 3.5” routine (Schulz and Mudelsee 2002). The appearance of the 1470-yr peak in the spectrum from the GISP2 data (Fig. 1) is caused by the regular spacing of three consecutive DO events (numbered 5, 6, and 7; Schulz 2002a). This regular spacing is not nearly as pronounced in the GRIP ice core using the flow model-based ss09–sea time scale (Johnsen et al. 1997), explaining why the spectral peak is not as significant in the GRIP spectrum. The discrepancy will hopefully be resolved in the near future where a more precise dating of the NGRIP ice core based on annual layers counting far back in the glacial period is expected. Levels of significance are not included since these are strongly depending on the null-hypothesis noise spectrum assumed. The regular Ornstein–Uhlenbeck process usually assumed (Schulz and Mudelsee 2002) is not relevant, since the isotope record is far from this process. Here we shall not argue for either the GISP2 or GRIP dating, although we note that one cannot favor one in comparison to the other solely based on the strength of the 1470-yr spectral peak as was done by Rahmstorf (2003).
Besides the uncertainty in dating, some discrepancy exists in defining the warming events. The “canonical” DO events (Dansgaard et al. 1993; Fig. 2) are based on visual inspection. Alley et al. (2001) used a filtering procedure, and threshold levels that produced 43 warming events were adopted. Schulz (2002a) found that DO events 9, 15, and 16 fell below the chosen threshold level, while the event around 65 kyr before present (BP) between events 18 and 19 was included. Rahmstorf (2003) only analyzed events younger than 50 kyr BP. Here event 9 was omitted while it was argued for an event (“A”) prior to the Younger Dryas (YD), and the termination of YD was included. We denote his event A as “GI1c” in accordance with Walker et al. (1999). In the following, we will argue for a procedure of defining the warming events. The problem is twofold: A spectral filtering must be decided, and threshold values and a procedure for threshold crossings must be chosen.
a. The filtering
There is a general consensus that a multimillennial high-pass filter must be applied in order to eliminate the variations due to the orbital forcing. Whether this is a spectral filter or subtraction of a long-term running mean is not significant. It is with respect to this long-term mean that the anomaly is defined. The important problem is to decide for a low pass in the other end of the spectrum. Most analyses have been performed on the approximately 200-yr-averaged isotope data available from the web archives. However, the spectral power in the time-scale range shorter than 200 yr is substantial. The effect of smoothing is that the positive extreme of a short time-scale warming event will be reduced and possibly brought below a chosen threshold level. Thus analyzing a smooth signal introduces a bias toward omitting very short DO events.
The shortest time-scale fluctuations in the isotope records are dominated by glaciological noise (Ditlevsen et al. 2002) and local fluctuations. The low pass should thus be determined such that the noise is reduced to an insignificant level in comparison to the true climatic fluctuations. By comparing the high-resolution GRIP and NGRIP records we can estimate the noise level. Figure 3 shows the correlation between the GRIP and the NGRIP records in the 27–37-kyr BP period, where the NGRIP signal has been dated by matching DO events 3–7 to the GRIP–ss09–sea time scale. The correlation coefficient is plotted as a function of the low-pass filter on both signals. The correlation between the raw high-resolution signals is approximately 0.5, while for the 30-yr low passes the correlation has increased to more than 0.8. The two core drilling sites are located more than 300 km apart so we conclude that the 30-yr low-pass signal is representative for climate fluctuations and thus appropriate for the analysis.
b. The threshold crossing
The two distinct quasi-stable climate states are most clearly seen in the bimodal distribution of the dust–calcium signal (Ditlevsen 1999) but is also apparent in the δ18O signal where the warm interstadial states are “sawtooth shaped,” characterized by gradual decreases in the signal prior to jumping into the stadial state (Alley 1998; Schulz 2002b). Accordingly, the signal can be split into these two states and two transition states, cold to warm and warm to cold (Ditlevsen et al. 2002). The initiation of the warm state is defined as the first up-crossing of the high threshold following an up-crossing of the low threshold. Similarly the initiation of the cold state is defined from the first down-crossing of the lower level following a down-crossing of the upper level (Fig. 4). The terminations are defined as the last up- (down) crossings of the lower (higher) level prior to an up- (down) crossing of the higher (lower) level. The periods between terminations and initiations are the transition periods.
The reason for this definition is that if the signal fluctuates around the higher (lower) level in one of the climatic states, this does not lead to a jump unless this follows a period where the signal was below (above) the lower (higher) level. This definition is advantageous in comparison to the definition used by Schulz (2002a). Here the signal is defined to be within the warm state whenever the value of the anomaly signal is above the higher threshold. This has the consequence that a warm period can more easily be split into more periods, as perhaps happens for DO event 17 (Schulz 2002a). The definitions applied here clearly identify whether DO event 17 should be regarded as one or two events. In our analysis, it remains one DO event (GI17).
In the following, the lower threshold is chosen as −1.0 permil anomaly and the higher threshold as +1.5 permil anomaly from the 10-kyr. high-pass isotope signal. The asymmetry is because the climate system persists in the cold state longer than in the warm state. The result is shown in Fig. 2. The vertical full lines indicate the first up-crossings into the warm states. The result is quite robust in the sense that the jumping times are rather insensitive to the threshold values chosen.
3. The waiting time distribution
The transitions from the cold to the warm state are much more abrupt than the opposite transitions, so the waiting times between consecutive up-crossings are less sensitive to the choices of filter frequencies and thresholds and are thus used for the further analysis. The cumulated distribution of waiting times is plotted in Fig. 5 (left). From a visual inspection, the observed distribution seems to be sampled from a Poisson process, which has an exponential distribution. The mean waiting time tm is 2.8 kyr, which is defined by the cumulated exponential distribution P(t) = 1 − exp(−t/tm) shown as the straight line. The discrepancy with the earlier findings (Alley et al. 2001; Schulz 2002a; Rahmstorf 2003) lies mainly in the inclusion of the shortest warming events. Whether or not the doubtful event “GS18b” is included does not change the waiting time distribution significantly and does not alter the following statistical analysis.
4. An upper bound on the periodic component
Even though the distribution shown in Fig. 5 (left) disfavors the previously claimed periodic component, the record is relatively short and, as Alley et al. (2001) noted, could be a realization of different possible processes. The waiting time distribution for a stochastic resonance process has a step-like structure with the first big step around the period and with exponentially smaller steps for multiples of the period (Fig. 5, right). A way of quantifying the degree of periodicity is then to calculate the root-mean-square (rms) difference between the distribution and the best-fit exponential distribution. Assuming the existence of a periodic component with period of 1470 yr, we can then estimate the strength of the periodic forcing in comparison to the barrier for purely noise-induced transition. Within the framework of the stochastic resonance model, this amounts to investigating the strength of the resonance.
We have generated a series of data from a stochastic resonance model with the same mean waiting time as observed in the isotope signal but with different strength of the periodic component. For each realization of length similar to the length of the isotope signal, we calculate the root-mean-square difference from the exponential distribution. The stochastic resonance model is described by the nonautonomous Langevin equation,
The first term represents the drift with two stable states separated by a potential barrier. The drift is derived from a potential F = −dU/dT, with U(T) = T4 − a3T3 − a2T2 + a1T, where a1, a2, and a3 are constants (Cessi 1994; Ditlevsen 1999). The second term in Eq. (1) is the periodic component with period τ and amplitude A. The third term is a white noise forcing with intensity σ.
For |A| ≥ Ac, Ac = −a1 − a2a3/2 − a33/8 + 8(a2/6 − a23/16)3/2, the system will jump between the two stable states through a hysteresis loop periodically, while for |A| < Ac the bifurcation points are not reached and the jumping must be noise induced.
A schematic of the potential and the periodic forcing is shown in Fig. 6. Disregarding nonexponential prefactors, the Kramer waiting time for penetrating the barrier going from the well containing the stable state a (c) to the well containing the stable state c (a) is estimated as Tac ∼ exp(Hmin/σ2) [Tca ∼ exp(Hmax/σ2)], where Hmin (Hmax) is the height of the barrier (Gardiner 1985). The criterion for the noise intensity to obtain stochastic resonance is
The climate state will then with high probability jump from state a to state c; within time τ/2 the potential has changed because of the periodic component, and the state will with high probability jump from c to a, and so on.
By varying the strength of the periodic component, and still tuning the noise intensity to the stochastic resonance, we can obtain an estimate of the amplitude of the periodic forcing component in comparison to the barrier in the system. Figure 7 shows a set of realizations of Eq. (1) for different strengths of the periodic forcing. Each realization is represented in a set of three panels composed of a long panel over two smaller panels. The strength of the stochastic resonance in the system is determined by the separation in time scales expressed in Eq. (2). When the noise is tuned to the resonance, the criterion is Δ ≡ Tca/Tac ∼ exp[(Hmax − Hmin)/σ2] ≫ 1. The four sets of panels correspond to Δ = 1.0, 1.6, 2.2, and 3.0, respectively.
Going from top to bottom, we see that the periodic component emerges. In the lower right panel of each set, the spectral peak at f = 2π/τ emerges and exceeds the 99% confidence level (Crowley et al. 1986) for Δ > 2.5. The cumulated waiting time distributions are shown in the left panels. The step-like structure from the periodicity emerges as Δ increases. The top set of panels (Δ = 1.0) shows the pure Poisson process. A quantitative measure of the deviation from the Poisson process is the rms distance of the waiting time distribution from an exponential. This is similar to the measure used for the Kolmogorov–Smirnov test, which cannot be applied when we obtain the best-fit distribution from the data. In the following, this rms is denoted the “error.” For each of the models shown in Fig. 7 the error is calculated from a sufficiently long simulation (Fig. 8). The accuracy is within the size of the plotting symbols.
To evaluate the data series, we have simulated a large ensemble of time series generated by Eq. (1) of same length as the isotope record for a given noise intensity. For each realization, we calculate the root-mean-square distance to the exponential distribution. By this procedure, the distribution of errors is calculated. In this way the maximum likelihood model for generating the observed record is obtained. Figure 9 shows the error distributions for the models generating the series shown in Fig. 7. The arrows represent the error obtained from the GRIP and GISP2 isotope records. The maximum likelihood model has Δ = 1 for the GRIP dating, while Δ = 1.6 for the GISP2 dating. For the model with Δ = 2.2, 2% of the realizations have an error less than the one measured for the isotope signal. We thus reject the hypothesis that the data are generated by a process with Δ ≥ 2.2 with 98% confidence for the GRIP dating. Correspondingly, we reject a process with Δ ≥ 3.0 with 98% confidence for the GISP2 dating. The NGRIP record is presently dated by fitting to the GRIP time scale and is thus not independent.
The isotope record shows the jumping between the warm DO interstadial state and the cold stadial state with waiting times in the millennial time-scale range. We have used an objective procedure based on correlating the GRIP and NGRIP records to decide the resolution of the climate record. A high and a low threshold for the anomalies are applied. The transitions are defined from consecutive level crossings. By this procedure we obtain the “canonical” DO events and an additional set of short events. This happens either by splitting some of the canonical events (1, 2, and 15) or by defining events previously in the cold periods (18 and 22). Assuming the record to be generated by a dynamics described by Eq. (1), the strength of a periodic component in the forcing in comparison to the strength of the barrier is expressed through the nondimensional parameter Δ. It will be relevant to extend the present analysis to include the oldest DO events from event 22 to the newly discovered event 25 (NGM) when a reliable dating is obtained for the NGRIP record. The observed record is highly probable as a realization of a purely noise-driven process without a periodic component. The shifts could be results of the erratic fluctuations in freshwater formation and heat transfer to the ocean surface. These fluctuations were strongest in the Last Glacial Maximum where there is a tendency for more frequent shifts (Rahmstorf 2003; Ditlevsen et al. 1996). The occurrence of event 25 shortly after the termination of the previous interglacial period (Eem) before substantial ice volumes have built up, suggests that the same climate modes exist in the interglacial periods where the “warm” mode persists (Ganopolski and Rahmstorf 2001). The reason why we do not experience DO events (except perhaps the 8.2-kyr event) in the Holocene climate could be the low intensity of fluctuations rather than the stability of the warm mode versus the “cold” mode as suggested by Ganopolski and Rahmstorf (2001).
Discussions with S. J. Johnsen are greatly appreciated. KKA thanks Carlsberg for funding.
Corresponding author address: Peter D. Ditlevsen, The Niels Bohr Institute, Dept. of Geophysics, University of Copenhagen, Juliane Maries Vej 30, DK-2100 Copenhagen O, Denmark. Email: firstname.lastname@example.org