Articles | Volume 19, issue 7
Research article
21 Jul 2023
Research article |  | 21 Jul 2023

Disparate energy sources for slow and fast Dansgaard–Oeschger cycles

Diederik Liebrand, Anouk T. M. de Bakker, Heather J. H. Johnstone, and Charlotte S. Miller

During the Late Pleistocene, Dansgaard–Oeschger (DO) cycles triggered warming events that were as abrupt as the present-day human-induced warming. However, in the absence of a periodic forcing operating on millennial timescales, the main energy sources of DO cycles remain debated. Here, we identify the energy sources of DO cycles by applying a bispectral analysis to the North Greenland Ice Core Project (NGRIP) oxygen isotope (δ18Oice) record; a 123 kyr long proxy record of air temperatures (Tair) over Greenland. For both modes of DO cyclicity – slow and fast – we detect disparate energy sources. Slow DO cycles, marked by multi-millennial periodicities in the 12.5 to 2.5 kyr bandwidth, receive energy from astronomical periodicities. Fast DO cycles, characterized by millennial periodicities in the 1.5 ± 0.5 kyr range, receive energy from centennial periodicities. We propose cryospheric and oceanic mechanisms that facilitate the transfer of energy from known sources to slow and fast DO cycles, respectively. Our findings stress the importance of understanding energy-transfer mechanisms across a broad range of timescales to explain the origins of climate cycles without primary periodic energy sources.

1 Introduction

Climate variability on (multi-) millennial timescales is an enigmatic phenomenon because there is no consensus about the energy sources for this variability (Stuiver and Braziunas, 1993; Wara et al., 2000; Schulz et al., 2002; Braun et al., 2005; Sun et al., 2021; Zhang et al., 2021). Astronomical climate forcing, i.e. a nonlinear response of the Earth's climate system to changes in the distribution of incoming solar radiation (insolation) across the hemispheres/seasons, can explain the variability in palaeoclimate time series in the Milankovitch band (i.e. tens to hundreds of thousands of years) (Hays et al., 1976; Liebrand and de Bakker, 2019; Riechers et al., 2022). Similarly, solar forcing through modest changes in the total annual amount of insolation the Earth receives, resulting from the  11-year sunspot cycle (±0.15 %), can explain palaeoclimate variability on timescales that range from the annual and decadal up to the centennial (Stuiver and Braziunas, 1993; Wagner et al., 2001; Peristykh and Damon, 2003; Huybers and Curry, 2006). The presence of energy sources for climate variability on these timescales starkly contrast with the absence of a widely accepted forcing agent on the in-between millennial and multi-millennial periodicities (i.e. thousands to tens of thousands of years) (Pestiaux et al., 1988; Hagelberg et al., 1991, 1994; Pelletier, 1998; Braun et al., 2005, 2009; Huybers and Curry, 2006; Scafetta et al., 2016; Rousseau et al., 2022). Insolation has been suggested to vary on these timescales (Scafetta et al., 2016; Kelsey, 2022), but the geological archive cannot be straightforwardly interpreted as a true recorder of (changes in the distribution of) solar output on millennial timescales due to the many nonlinear response mechanisms of the Earth's climate system. To advance the debate about what “fuels” millennial climate cycles, it is thus essential to see whether millennial climate variability is cyclic (i.e. periodic) (Ditlevsen et al., 2005), and if so, where the energy is derived from.

The most notable example of millennial climate variability is the Dansgaard–Oeschger (DO) cycles observed in the North Greenland Ice Core Project (NGRIP) oxygen isotope (δ18Oice) record (Andersen et al., 2004). As an air temperature (Tair) proxy record for the past 123 kyr (Late Pleistocene to Holocene), the NGRIP δ18Oice chronology is unsurpassed in its resolution, dating precision and accuracy, and detail of structure in the data (i.e. signal-to-noise ratio) (Andersen et al., 2004, 2006; Rasmussen et al., 2006; Vinther et al., 2006; Svensson et al., 2008; Wolff et al., 2010; Kindler et al., 2014). This level of data quality may constitute a permanent ceiling in what can be reconstructed using climate proxies from any of the Earth's natural archives, because no marine or terrestrial record will probably ever be able to match ice core chronologies in all these respects simultaneously. A comprehensive statistical analysis of the NGRIP δ18Oice record – including a higher-order spectral analysis – is thus essential if we want to fully probe the climate dynamics across centennial, millennial, and astronomical timescales captured in this unique chronology.

Previous statistical tests showed that certain spectral peaks in marine records were lower or higher harmonics of a primary cyclic signal (Hagelberg et al., 1994; Riechers et al., 2022), findings suggestive of frequency- and phase-coupling among climate cycles. However, these analyses were not able to detect direction and magnitude of energy transfers among climate cycles. We overcome the shortcomings of these methods by computing directions and magnitudes of energy transfers among climate cycles in the NGRIP δ18Oice record directly from the bispectrum – a statistical technique originally developed by the climatologist and 2021 Nobel Prize laureate in Physics Klaus Hasselmann and colleagues to investigate the nonlinear origins of near-shore ocean waves (Hasselmann et al., 1963). Our interpretation of the bispectral results is based on the premise that the Earth's climate–cryosphere system behaves, at least partially, as a nonlinear oscillator (Le Treut and Ghil, 1983). We apply an advanced interpretation of the bispectrum to the NGRIP δ18Oice record (i.e. the total and zonal integrations of the imaginary part of the bispectrum) (Herbers et al., 2000; de Bakker et al., 2015; Liebrand and de Bakker, 2019) and show that both slow and fast DO cycles derive their energy from largely disparate astronomical and centennial sources through nonlinear interactions inherent to the climate–cryosphere system. These findings may abolish the need for invoking millennial forcing agents (Scafetta et al., 2016) to explain millennial climate cycles.

2 Asymmetric DO cycles in the NGRIP δ18Oice record

In total, 26 quasi-periodic DO cycles are currently recognized in the NGRIP δ18Oice record (see Sect. 6), some of which are subdivided into a few shorter lasting episodes (Fig. 1a) (Dansgaard et al., 1993; Andersen et al., 2004; Rasmussen et al., 2014). They are characterized by abrupt warming events at the base of relatively warm interstadials (δ18Oice=−36 ‰ to −40 ‰, Tair-35 to −40C), followed by either several stepwise and rapid, or one more gradual, Tair decrease(s), which lead(s) to longer-lasting and colder stadials (δ18Oice=-40 ‰ to −45 ‰, Tair-50C) (Fig. 1a) (Kindler et al., 2014). For most DO events, this evolution, i.e. an abrupt warming followed by a stepwise or a more gradual cooling, results in asymmetric (i.e. sawtooth-shaped) cycle shapes or waveforms (Hagelberg et al., 1991; King, 1996). According to bispectral theory, asymmetric cycle shapes in time series constitute a prognostic geometry for the nonlinear transfer of energy across the power spectrum (Herbers et al., 2000; de Bakker et al., 2015).

Outside of Greenland, climate variability that is highly comparable to the DO cycles in structure is recognized in planktonic foraminiferal δ18O values from the Iberian margin, indicating that sea surface temperatures in the wider North Atlantic region were as variable on millennial timescales as the Tair cycles over Greenland (Shackleton et al., 2000; Hodell et al., 2023). Evidence for the global climatic expression of Late Pleistocene millennial climate variability, albeit more subdued, comes from, e.g. benthic foraminiferal δ18O values (i.e. deep water temperatures) from the Iberian margin (Shackleton et al., 2000; Hodell et al., 2023), δD (i.e. air temperature) and methane records from the Antarctic ice core chronology (Loulergue et al., 2008), and globally distributed speleothem records (Corrick et al., 2020; Menviel et al., 2020; Batchelor et al., 2023).

Figure 1The NGRIP δ18Oice record versus summer insolation at 65 N. (a) The NGRIP δ18Oice record on the “GICC05modelext” age model (Andersen et al., 2006; Rasmussen et al., 2006; Vinther et al., 2006; Svensson et al., 2008; Wolff et al., 2010) versus (b) mean summer half-year insolation at 65 N (21 March to 21 September) for the past 123 kyr (Laskar et al., 2004). (a) The grey line represents the NGRIP δ18Oice record with 20-year sample spacing, used for bispectral analyses (Figs. 2–4). The black line represents the NGRIP δ18Oice record with a 50-year sample spacing, used for the wavelet analysis. (c) Wavelet analysis of the Gaussian filtered NGRIP δ18Oice record (see Sect. 6). Black contours represent the 95 % significance level with respect to a red noise model. The colours of the wavelet indicate amplitude (white/grey = low amplitude, yellow/orange/red = high amplitude). We note that colour/amplitude does not equate to statistical significance with respect to the red noise model. An average of the entire wavelet spectrum is shown on the right-hand side in red. The two dashed white lines indicate the spectral gaps between astronomical periodicities and slow DO periodicities (broad dashes), and between slow DO periodicities and fast DO and centennial periodicities (narrow dashes); “H” stands for Holocene, “BA” stands for Bølling–Allerød, and “YD” stands for Younger Dryas. The latter two together correspond to DO cycle 1.


Prior to investigating the asymmetric cycle shapes of the NGRIP δ18Oice record in more detail, we reassess the (quasi-) periodic nature of millennial climate cycles (Fig. 1c) (Schulz et al., 1999; Alley et al., 2001; Hinnov et al., 2002; Matyasovszky, 2010). We note that there is a lack of consensus in the literature about the presence/absence of periodic components in the NGRIP δ18Oice record, with some studies arguing against any form of periodicity but in favour of purely stochastically induced DO “events” (Ditlevsen et al., 2005, 2007; Kleppin et al., 2015; Gottwald, 2021). However, considering our new bispectral results that identify a clear grouping in the distribution of nonsinusoidality (i.e. asymmetry) per frequency (see next paragraph), we also believe that a reassessment of the distribution of sinusoidality per periodicity is justified. To do so, we apply a wavelet analysis to the data (Fig. 1c) (see Sect. 6). A Gaussian notch filter, applied to the NGRIP δ18Oice record prior to wavelet analysis, somewhat enhances the power of the millennial periodicities with respect to the astronomical periodicities. Similar to previous studies that performed spectral analysis on the NGRIP δ18Oice record (Petersen et al., 2013; Mitsui et al., 2019; Riechers et al., 2022), we identify two distinct modes of DO variability in the time–periodicity domain: (i) slow DO cycles, which are marked by sub-astronomical/multi-millennial periodicities in the 12.5 to 2.5 kyr bandwidth (e.g. DO cycles 1, 8, 12, 19, 20, and 22) (Dansgaard et al., 1984) and (ii) fast DO cycles that have millennial periodicities in the 1.5 ± 0.5 kyr range (Dansgaard et al., 1984; Schulz et al., 1999; Schulz, 2002; Braun et al., 2005; Ditlevsen et al., 2007) that extend into the centennial periodicity range (e.g. DO cycles 3, 4, 5, 6, 7, 9, 10, 11, and 18) (Fig. 1c). Slow DO cycles are near-continuously present throughout the record and are especially strongly expressed from 110 to 105 ka, from 80 to 70 ka, from 50 to 35 ka, and from 15 to 10 ka. Slow DO cycles contribute to the variance associated with the Eemian and Holocene interglacials, as well as high amplitude DO cycles 1 (i.e. Bølling–Allerød/Younger Dryas), 19, and 20 (Fig. 1). Fast DO cycles are briefly present from 110 to 105 ka, and near-continuously expressed from 90 to 10 ka. The distinction in cycle durations between the two modes of DO variability suggests that there may be at least two distinct climatic mechanisms at their root cause.

3 Energy transfers among climate cycles shed light on the “black box” climate system

To further probe the origins of both slow and fast DO cycles in the NGRIP δ18Oice record, we apply a bispectral analysis (Fig. 2) to the detrended and Hamming-tapered data (see Sect. 6). The Hamming tapering focusses the time-averaged bispectral analysis on the central part of the 123 kyr long NGRIP δ18Oice record. Bispectral analysis is a statistical technique that can deconvolve the asymmetry (i.e. sawtooth shape) of a periodic signal into three parts. These constituents, one sum frequency (f3, i.e. f3=f1+f2) and two difference frequencies (f1 and f2, i.e. f1=f3-f2, and f2=f3-f1), can either gain (orange and red colours for f3, blue colours for f1 and f2) or lose (blue colours for f3, orange and red for f1 and f2) energy (Figs. 2, S1, S2, Table S1). This makes bispectral analysis a powerful technique to identify energy sources of asymmetric periodic signals. In the bispectrum of the NGRIP δ18Oice record (Fig. 2), we observe three slow DO periodicities that stand out by being involved in many energy exchanges, namely the 6.1, 4.2, and 3.7 kyr cycles. For fast DO cycles, we observe more subtle interactions at the crossing lines with the 40 kyr obliquity and 23 kyr precession cycles, and, more prominently, with the 4.2 and 3.7 kyr slow DO cycles. Relative energy transfers among climate cycles, and their directions, can thus be derived from so-called triad interactions displayed in (the imaginary part of) the bispectrum. The bispectrum shows that slow and fast DO cycles have highly complex and mixed energy sources from across a wide range of periodicities. However, the net effect per periodicity of all individual energy exchanges cannot be derived from the bispectrum alone because individual periodicities are represented in the bispectrum as f1, f2, and f3 simultaneously (i.e. 1/p1, 1/p2, and 1/p3), each of which can gain and/or lose energy (see Sect. 6) (de Bakker et al., 2015; Liebrand and de Bakker, 2019).

Figure 2Bispectrum of the NGRIP δ18Oice record. The imaginary part of the bispectrum of the NGRIP δ18Oice record depicts energy exchanges among three frequencies. Energy can only be transferred if both frequencies (f) and phases (φ) are coupled: f1+f2=f3 and φ1+φ2=φ3. The x axis shows difference frequency f1, the y axis depicts difference frequency f2, and the colour bar corresponds to sum frequency f3 (i.e. f1+f2). Two outcomes exist for energy exchanges: either f3 gains energy and f1 and f2 lose energy, which is shown in orange and red colours, or f3 loses energy and f1 and f2 gain energy, which is shown in blue colours. Purple-, plum-, and pink-coloured bands mark the main periodicities for astronomical, slow, and fast DO cycles, respectively. Vertical and horizontal coloured bands mark f1 and f2, respectively. Diagonal bands (curved on the log–log scale) mark f3. Thick solid black lines mark our selected zone boundaries between astronomical periodicities, slow DO periodicities, and fast DO and centennial periodicities. See Fig. S2 in the Supplement for a linear–linear scale version of this figure. For a more detailed explanation of bispectral analysis, we refer to Sect. 6 and Liebrand and de Bakker (2019).


To obtain magnitudes of net energy transfers per periodicity (i.e. “climate cycle”), we need to integrate over the bispectrum (see Sect. 6). We apply a coupling coefficient as part of this integration step to scale energy with frequency. The total integration of the bispectrum of the NGRIP δ18Oice record reveals three unique bandwidths in the periodicity domain, each marked by distinct regimes in energy losses and/or gains (Fig. 3b). First, the astronomical periodicities of ∼110 kyr eccentricity, 40 kyr obliquity, and ∼20 kyr precession, lose energy. Second, periodicities of slow DO cycles, covering the bandwidth from 12.5 to 2.5 kyr, predominantly gain energy. Third, fast DO cycles and centennial periodicities, in the 2.5 to 0.1 kyr periodicity range, are marked by rapidly alternating net energy gains and losses that seem to largely balance one another across this part of the spectrum, notwithstanding an overall loss in energy at periodicities in the range between 0.3 and 0.1 kyr (Fig. 3b). The net loss of energy at the ∼110 kyr eccentricity periodicity is a surprising result that contrasts with a Milankovitch-Theory-based understanding of precession and obliquity fuelling of eccentricity variability of the Earth's climate–cryosphere system (Liebrand and de Bakker, 2019); this is a result that should be interpreted with caution (see Sect. 6). Net energy gains that stand out in amplitude are the slow DO cycles with periodicities of 6.1, 4.2, and 3.7 kyr, as well as the gains at the 2.9 and 2.7 kyr periodicities. Some further notable net energy gains coincide with the 1.5 ± 0.5 kyr periodicity range of the fast DO cycles, such as the 1.7, 1.5, and 1.3 kyr periodicities. Other well-expressed gains fall in the centennial bandwidth, being the 0.9, 0.5, 0.4, and 0.2 kyr periodicities. By integrating over the bispectrum, we show which periodicities gain and lose energy, and thus how energy is redistributed over the power spectrum. Based on this new information, we can piece together cause-and-effect relationships otherwise unobservable within the “black box” of the Earth's climate system and derive an origin story for both slow and fast DO cycles.

Figure 3Linking DO cycles to energy sources. Comparison of (a) power spectrum of mean summer half-year insolation record at 65 N (21 March to 21 September) to (b) total integration of the bispectrum of the NGRIP δ18Oice record and (c) the power spectrum of the NGRIP δ18Oice record. Orange indicates spectral power (a, c), or net energy gains (b), and blue indicates net energy losses (b). Red and orange arrows indicate the two main energy-transfer pathways across the spectral continuum to slow and fast DO cycles. These interpretative arrows are derived from both the total (this figure) and the zonal integration (see also Fig. 4). Dashed vertical black lines mark the zone boundaries between astronomical periodicities (purple), slow DO cycles (plum), and fast DO cycles (pink) and centennial periodicities; “S” stands for sun/solar forcing, indicating the frequencies at which summer insolation provides energy. There is no 123 kyr long record of sunspot cycles, which is why we do not show a spectrum of this energy source here. Furthermore, the highest amplitude spectral power of the solar cycle over the past 123 kyr is expected to fall at around the annual and ∼11-year periodicities. These cycles are not resolved in the NGRIP δ18Oice record, which is marked by a 20-year sample spacing.


By comparing the power spectrum of insolation (Fig. 3a) to the total integration of the bispectrum (Fig. 3b) and power spectrum (Fig. 3c) of the NGRIP δ18Oice record, we show that part of the energy the Earth receives on astronomical timescales (Fig. 3a) is lost (Fig. 3b) and transferred toward shorter periodicities, particularly those of slow DO cycles (Fig. 3b, c). Slow DO cycles thus have an indirect astronomical origin (Hagelberg et al., 1994; Wara et al., 2000; Rial and Saha, 2011; Sun et al., 2021; Zhang et al., 2021), at least in part. Furthermore, we observe several distinct energy gains in the 1.5 ± 0.5 kyr periodicity range of fast DO cycles, which are accompanied by smaller amplitude energy losses. These net gains show that part of the spectral power associated with fast DO cycles is derived from nonlinear interactions with other climate cycles (Fig. 3b, c). However, from the total integration of the bispectrum, it is hard to make out whether most of this energy is derived from astronomical or other fast DO and centennial cycles.

To better identify the disparate energy sources for both slow and fast DO cycles, we devise a zonation scheme (Fig. 2). This zonation scheme aims to isolate groups of nonlinear interactions among specific frequency bandwidths. It consists of nine unique combinations of the three main bandwidths: (i) astronomical, (ii) slow DO, and (iii) fast DO and centennial periodicities (Figs. 2, S1, Table S1). Boundaries for the zonal integration are selected at frequencies where the qualitative behaviour of energy transfers in the total integration changes (Fig. 3b). This results in f=80 Myr−1 (p=12.5 kyr) for the boundary between astronomical and slow DO cycles, and at f=400 Myr−1 (p=2.5 kyr) for the boundary between slow DO and fast DO (and centennial) cycles. Subsequently, we integrate the bispectrum of the NGRIP δ18Oice record zonally (Fig. 4) (see Sect. 6) (de Bakker et al., 2015; Liebrand and de Bakker, 2019). From the zonal integrations, we learn that most energy is exchanged among astronomical and slow DO cycles (Zone 3, Fig. 4c), which results in a net transfer of energy to the 3.7 kyr cycle. Other energy sources for slow DO cycles now also become apparent; additional energy is sourced from fast DO periodicities in the range between 2.0 and 2.5 kyr (Zone 6, Fig. 4f) and broadly from across the millennial to centennial periodicity range (Zone 8, Fig. 4h). Fast DO cycles, on the contrary, do not appear to receive much energy from astronomical periodicities (Zone 7, Fig. 4g) but instead are “fuelled” by a combination of centennial climate cycles and secondary interactions with other fast DO periodicities (Zones 6–9, Fig. 4f–i).

Figure 4Zonal integrations of the bispectrum. Integrations of (a) Zone 1: BIm(A, A, A), (b) Zone 2: BIm(A, A, DOS), (c) Zone 3: BIm(DOS, A, DOS), (d) Zone 4: BIm(DOS, A, DOF), (e) Zone 5: BIm(DOS, DOS, DOS), (f) Zone 6: BIm(DOS, DOS, DOF), (g) Zone 7: BIm(DOF, A, DOF), (h) Zone 8: BIm(DOF, DOS, DOF), and (i) Zone 9: BIm(DOF, DOF, DOF); “BIm” stands for the imaginary part of the bispectrum, “A” stands for astronomical periodicities, “DOS” stands for slow DO periodicities, “DOF” stands for fast DO and centennial periodicities. Orange indicates net energy gains while blue represents net energy losses. Dashed vertical black lines mark the zone boundaries between the main periodicity bandwidths, with the main astronomical periodicities highlighted in purple, slow DO periodicities in plum, and fast DO periodicities in pink.


4 Potential cryospheric, climatic, and oceanic energy-transfer mechanisms

By statistically reanalysing the NGRIP δ18Oice record, we find clear confirmatory (Petersen et al., 2013; Vettoretti and Peltier, 2018) evidence for two distinct modes (i.e. slow and fast) of DO variability (Figs. 1c and 3b). Furthermore, based on our advanced bispectral analyses (Figs. 2, 3b, and 4), we show that disparate fuelling pathways underpin these two modes of DO cyclicity. If these interpretations are correct, then slow DO cycles constitute a nonlinear response primarily to astronomical climate forcing (Hagelberg et al., 1994; Wara et al., 2000; Rial and Saha, 2011; Sun et al., 2021; Zhang et al., 2021), which, in combination with secondary energy-redistributing mechanisms among slow DO periodicities, results in Tair cyclicity over Greenland in the 12.5 to 2.5 kyr bandwidth. It is important to note that the direction of energy transfers, from obliquity and precession periodicities (and potentially eccentricity periodicities; see Sect. 6) to the multi-millennial slow DO periodicities (Figs. 2, 3b, and 4), is inversed (i.e. from longer to shorter periodicities) from what has been documented for the Middle and Late Pleistocene 40 kyr and ∼110 kyr worlds (Liebrand and de Bakker, 2019). These ice age cycles are, in accordance with the Milankovitch Theory, marked by energy transfers from precession to obliquity periodicities and from precession and obliquity to eccentricity periodicities, respectively (i.e. from shorter to longer periodicities). We speculate that the energy-redistributing mechanisms from precession and obliquity periodicities to and among slow DO periodicities were predominantly linked to the Northern Hemisphere (NH) climate–cryosphere system (Vettoretti and Peltier, 2018). For example the regeneration time, or waiting time (Alley et al., 2001), of grounded ice in the Hudson Bay following catastrophic collapses (Boers et al., 2018) could be astronomically forced. In such a scenario, the rate of slow DO cooling phases and grounded ice expansion (i.e. inceptions), potentially in combination with floating ice shelf expansions near Greenland, were paced by insolation conditions. Yet their warming phases (i.e. terminations), would have been triggered by the passing of (an) internal threshold(s), most probably during the NH summer. Such a mechanism, somewhat comparable to the binge/purge model for continental-sized ice sheets (MacAyeal, 1993), would result in asymmetric cycle shapes that mark the slow DO cycles, as well as have a short periodic nonlinear response to insolation forcing built in, as we document in our bispectral results. These quasi-regular oscillations would have affected sea-ice cover and regional Tair, most probably also over Greenland. Furthermore, the relatively long time periods needed for ice sheets to respond are in line with a relatively slow mode of DO cyclicity.

For fast DO cycles, the energy redistributing mechanisms (i.e. from shorter to longer periodicities) are not as straightforward as those for slow DO cycles. However, our results suggest that many, closely spaced energy exchanges in the frequency domain (see e.g. the rapidly alternating losses and gains in Fig. 3b) may have served as a primary energy source for fast DO cycles in the 1.5 ± 0.5 kyr bandwidth. Although energy for fast DO cycles may have been collected from many centennial periodicities (Ditlevsen and Crucifix, 2017), mainly in the 0.3 to 0.1 kyr range (Figs. 3b and 4), initially this energy was probably derived from long-term modulations in the solar cycles (Wagner et al., 2001; Peristykh and Damon, 2003), or derived from even shorter periodic (e.g. decadal) sunspot cycles not captured in our analysis. We tentatively support a solar forcing hypothesis for fast DO cycles, albeit not through a summing of periodicities (Braun et al., 2005), but through nonlinear energy exchanges (i.e. triad interactions) with shorter periodicities. We speculate that changes in ocean currents linked to Atlantic meridional overturning circulation (AMOC) and/or tides may be the primary driver for the distinct fast DO mode of variability (van Kreveld et al., 2000; Elliot et al., 2002; Schmidt et al., 2006; Dokken et al., 2013; Vettoretti and Peltier, 2018; Griem et al., 2019; Armstrong et al., 2022). In this scenario, AMOC strength had destabilizing effects on floating ice shelves only, but it did not affect grounded ice as much. Ice shelf collapse and rapid warming followed by a somewhat slower phase of ice shelf re-expansion would explain the asymmetric shape of the fast DO cycles (Schulz et al., 1999). Both modes of DO variability, slow and fast, were likely amplified through positive feedbacks of the carbon cycle (Brook et al., 1996; Bauska et al., 2021; Vettoretti et al., 2022). In fact, by chemically storing energy at certain timescales/periodicities, and releasing it at others, the carbon cycle itself could have functioned as an energy-redistributing mechanism. A rechargeable battery with different charging and discharging rates may serve as a metaphor for such a process. It is thus important to consider both physical (e.g. build up and loss of grounded and floating ice) and chemical (e.g. energy stored in, and released from, organic carbon) energy sinks and sources for the nonlinear redistribution of energy across timescales, as quantified in the bispectrum.

5 Conclusions

DO cycles, both slow and fast, resulted – at least in part – from cryospheric, climatic, and oceanic energy-redistributing processes that operated on a broad continuum of periodicities and linked dynamic behaviour across timescales (Schulz et al., 1999, 2002; Schulz, 2002; Markle et al., 2017; Boers et al., 2018; Menviel et al., 2020; Armstrong et al., 2022; Vettoretti et al., 2022). We provide a description of two distinct energy-transfer pathways: one from longer astronomical periodicities and one from shorter centennial periodicities to the (multi-) millennial periodicities that mark the DO cycles. We emphasize that the disparate energy sources for DO variability that we identify in the bispectrum of the NGRIP δ18Oice record constitute periodic Earth internal sources of energy. However, both these Earth internal sources are in turn externally forced by insolation variability on either astronomical or centennial (and shorter) timescales. The disparate energy sources and associated climate mechanisms contrast with previous studies that favour a single mechanism to explain both slow and fast DO cyclicity (e.g. Petersen et al., 2013; Lohmann and Svensson, 2022; Vettoretti et al., 2022). The bispectrum can only assess net energy exchanges among periodic components of a time series, and hence, we cannot rule out non-periodic, non-energy-redistributing forcing factors, such as volcanism (Crick et al., 2021; Paine et al., 2021; Lin et al., 2022; Lohmann and Svensson, 2022), asteroid impacts (van Hoesel et al., 2014), and/or noise-induced or other “stochastic events” (Alley et al., 2001; Braun et al., 2009; Lohmann et al., 2020; Vettoretti et al., 2022). However, our study stresses the importance of considering nonlinear – but deterministic – energy-redistributing mechanisms within the climate–cryosphere system to explain periodic signals on timescales without primary periodic energy sources.

6 Methods

NGRIP age model. The NGRIP record is composed of two bore holes that are spliced together. These cores recovered snow and ice layers with annual resolution. To generate a high-fidelity time series, 20 layers were averaged, yielding a highly reproducible and largely noise-free δ18Oice chronology. The ice-water δ18O record constitutes a Tair proxy record of the past 123 kyr that captures astronomical, millennial, and centennial climate variability in exquisite detail (Andersen et al., 2004; Kindler et al., 2014) (Fig. 1). The “GICC05modelext” age model (Andersen et al., 2006; Rasmussen et al., 2006; Vinther et al., 2006; Svensson et al., 2008; Wolff et al., 2010) used here was constructed using a combination of annual ice layer counting and an ice flow model. We note that the relatively small age model differences (up to ∼0.5 kyr) compared to the older “ss09sea” age model (Andersen et al., 2004) do not affect the outcome of bispectral analyses because bispectra assess cycle geometry in the frequency domain, which allows for modest changes in the time domain (Liebrand and de Bakker, 2019).

Statistical analyses. The wavelet analysis of the NGRIP δ18Oice record transforms the time series into the time–periodicity domain using a Morlet wavelet (Fig. 1c). Prior to the wavelet analysis only, we applied a Gaussian notch filter (f=0 kyr−1, bandwidth = 0.05 kyr−1, Paillard et al., 1996) to the NGRIP δ18Oice record (50-year sample spacing) to somewhat enhance the power of the millennial periodicities with respect to the astronomical periodicities (i.e. the latter are suppressed by the notch filter).

In this study, we use the same bispectral methods as described in Liebrand and de Bakker (2019). In short, the bispectrum is defined by Eq. (1):

(1) B f 1 , f 2 = E [ C f 1 C f 2 C f 1 + f 2 * ] ,

where E[] is the ensemble average of the triple product of complex Fourier coefficients C at the difference frequencies f1, f2, and their sumf3 (i.e. =f1+f2), and the asterisk indicates complex conjugation (Hasselmann et al., 1963). Energy gains and losses in both the bispectrum and the total and zonal integration of the bispectrum are inverted compared to Liebrand and de Bakker (2019) because δ18Oice records an inverted ice age signal compared to benthic foraminiferal δ18O. To extract the total gain or loss in energy per climate cycle over the 123 kyr span of the NGRIP δ18Oice record, and hence, to obtain conservative net energy transfers per frequency (defined by the nonlinear source term Snl), we integrate over the imaginary part of the bispectrum and multiply it with a coupling coefficient. Similar to Liebrand and de Bakker (2019), we make minimum assumptions and use a coupling coefficient that only corrects for frequency W(f1,f2)=(f1+f2), which is also part of the Boussinesq scaling used for ocean waves (Herbers and Burton, 1997; Herbers et al., 2000). The correction for frequency (W(f1,f2)) ensures energy conservation during triad interactions (see Fig. S3) and permits qualitative interpretations of conservative net energy exchanges across the spectrum (see Figs. 3b and 4). We express the integral of Snl in terms of an integration, or summation, over the positive quadrant of the bispectrum alone, or equivalently in sum and difference interactions following Eq. (2):

(2) S nl , f = S nl , f + + S nl , f - .

The difference contributions are expressed by Eq. (3):

(3) S nl - = - 2 f F - W f + f , - f I { B f , f } ,

and are obtained by integrating along vertical and horizontal lines, perpendicular to the x and y axes, respectively (Figs. 2, S2). The sum contributions are expressed by Eq. (4):

(4) S nl + = f F + W f , f - f I { B f , f - f } ,

and are obtained by integrating diagonally over the bispectrum (Figs. 2, S2). Integration can be performed including all frequencies, i.e. a total integration over the entire imaginary part of the bispectrum (Fig. 3b), or over specific zones only, in case subsets of interactions need to be further examined (Fig. 4, Table S1) (de Bakker et al., 2015, 2016).

The imaginary part of the bispectrum, the total, and the zonal integrations constitute an average of the past 123 kyr. Prior to bispectral analysis, the raw NGRIP δ18Oice record (20-year sample spacing) was linearly detrended and multiplied by a Hamming window function. We subsequently used an energy correction factor that corrects for the Hamming tapering. Previous statistical analyses (Liebrand and de Bakker, 2019) have shown that this approach considerably clarifies the bispectral results. Resulting from the tapering, the central part of the time series is relatively overrepresented, whereas those parts on either end of the record (i.e. the Eemian, Bølling–Allerød/Younger Dryas, and Holocene) are relatively underrepresented in the bispectral analysis. Furthermore, only one ∼110 kyr cycle is present in the 123 kyr long NGRIP δ18Oice record. The structure of this ∼110 kyr cycle, and hence the nonlinear triad interactions in the bispectrum, are affected by the Hamming tapering. We thus caution against overinterpreting the energy loss observed at the ∼110 kyr eccentricity periodicity (i.e. the very lowest frequency), which is at, or beyond, the detection limit of our analysis. Future analyses on longer records, marked by both astronomical and millennial climate variability, are needed to verify the robustness of the energy loss at the ∼110 kyr eccentricity periodicity in the NGRIP δ18Oice record, because longer window lengths, that incorporate multiple ∼110 kyr cycles, can then be selected (see also Liebrand and de Bakker, 2019). We refer the interested reader to Schmidt (2020) for some of the latest applications and interpretations of bispectral analysis.

Following studies on ocean waves (Herbers et al., 2000; de Bakker et al., 2015), we focus the interpretation on the imaginary part of the bispectrum (Figs. 2, 3b, and 4), which is related to cycle asymmetry, because most energy transfers are associated with asymmetric cycle shapes/wave forms. Furthermore, we visually observe sawtooth-shaped/asymmetric climate cycles in the NGRIP δ18Oice record, associated with both slow and fast DO cycles. It is worth noting that the asymmetry/sawtooth shape of DO cycles, analysed here using the imaginary part of the bispectrum, is distinct from their skewness/peakedness (Ditlevsen et al., 2005), which is deconvolved in the real part of the bispectrum, or from excess kurtosis/rectangularity/square-waveness (Ditlevsen et al., 1996; Hinnov et al., 2002), which is deconvolved in the trispectrum. The latter two constitute wholly different nonsinusoidal cycle properties (Hagelberg et al., 1991; King, 1996). Skewed and kurtose cycle geometries are not as straightforwardly associable with energy transfers (Herbers et al., 2000; de Bakker et al., 2015; Liebrand and de Bakker, 2019).

The bispectrum yields nonconservative relative energy transfers that have not yet been corrected for frequency (de Bakker et al., 2015; Liebrand and de Bakker, 2019). Prior to such a scaling, the relative energy transfers at higher frequencies are not as pronounced as those at lower frequencies in the bispectrum (Fig. 2). During the integration of the bispectrum, we scale energy transfers to frequency by assuming a standard coupling coefficient, which is based on the second order Boussinesq theory (Herbers and Burton, 1997; Herbers et al., 2000). This scaling yields conservative net energy transfers for the imaginary part of the bispectrum only (de Bakker et al., 2015; Liebrand and de Bakker, 2019) (Figs. 3, 4, and S3). We note, however, that both the nonconservative relative energy transfers of the bispectrum, as well as the conservative net energy transfers of the integration of the bispectrum, are not absolute (Liebrand and de Bakker, 2019). We forego scaling energy transfers to insolation due to the many unknowns in the climate–cryosphere system that link changes in δ18Oice (i.e. in ‰) to changes in insolation (i.e. in W m−2).

Code availability

The detailed mathematical description of the bispectral analyses used in this study are presented in Sect. 6 as well as in de Bakker et al. (2015, and Liebrand and de Bakker (2019, Code will be made available on individual request to the corresponding authors.

Data availability

The NGRIP δ18Oice record can be downloaded here: (University of Copenhagen, Centre for Ice and Climate, Niels Bohr Institute, 2023).


The supplement related to this article is available online at:

Author contributions

DL and ATMdB designed the study, performed the data analysis, and made the figures. ATMdB adjusted the bispectral MATLAB scripts for the purposes of this study. All authors discussed the interpretation of the results and contributed to the writing of the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “Ice core science at the three poles (CP/TC inter-journal SI)”. It is not associated with a conference.


We thank the NGRIP members and the wider ice-core drilling and science community for making the δ18Oice data and accompanying age models available. We thank Thomas H. C. Herbers and B. Gerben Ruessink for making their original bispectral MATLAB scripts available, which were adapted for this study. We thank Andrew Yool for help with the three-dimensional visualization of the wavelet analysis.

Financial support

The article processing charges for this open-access publication were covered by the University of Bremen.

Review statement

This paper was edited by Alexey Ekaykin and reviewed by David De Vleeschouwer and one anonymous referee.


Alley, R. B., Anandakrishnan, S., and Jung, P.: Stochastic resonance in the North Atlantic, Paleoceanography, 16, 190–198,, 2001. 

Andersen, K. K., Azuma, N., Barnola, J. M., Bigler, M., Biscaye, P., Caillon, N., Chappellaz, J., Clausen, H. B., Dahl-Jensen, D., Fischer, H., Fluckiger, J., Fritzsche, D., Fujii, Y., Goto-Azuma, K., Gronvold, K., Gundestrup, N. S., Hansson, M., Huber, C., Hvidberg, C. S., Johnsen, S. J., Jonsell, U., Jouzel, J., Kipfstuhl, S., Landais, A., Leuenberger, M., Lorrain, R., Masson-Delmotte, V., Miller, H., Motoyama, H., Narita, H., Popp, T., Rasmussen, S. O., Raynaud, D., Rothlisberger, R., Ruth, U., Samyn, D., Schwander, J., Shoji, H., Siggard-Andersen, M. L., Steffensen, J. P., Stocker, T., Sveinbjornsdottir, A. E., Svensson, A., Takata, M., Tison, J. L., Thorsteinsson, T., Watanabe, O., Wilhelms, F., White, J. W. C., and Project, N. G. I. C.: High-resolution record of Northern Hemisphere climate extending into the last interglacial period, Nature, 431, 147–151,, 2004. 

Andersen, K. K., Svensson, A., Johnsen, S. J., Rasmussen, S. O., Bigler, M., Rothlisberger, R., Ruth, U., Siggaard-Andersen, M. L., Steffensen, J. P., Dahl-Jensen, D., Vinther, B. M., and Clausen, H. B.: The Greenland Ice Core Chronology 2005, 15–42 ka. Part 1: constructing the time scale, Quaternary Sci. Rev., 25, 3246–3257,, 2006. 

Armstrong, E., Izumi, K., and Valdes, P.: Identifying the mechanisms of DO-scale oscillations in a GCM: a salt oscillator triggered by the Laurentide ice sheet, Clim. Dynam., 60, 3983–4001,, 2022. 

Batchelor, C. J., Marcott, S. A., Orland, I. J., He, F., and Edwards, R. L.: Decadal warming events extended into central North America during the last glacial period, Nat. Geosci., 16, 257–261,, 2023. 

Bauska, T. K., Marcott, S. A., and Brook, E. J.: Abrupt changes in the global carbon cycle during the last glacial period, Nat. Geosci., 14, 91–96,, 2021. 

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, 47,, 2018. 

Braun, H., Christl, M., Rahmstorf, S., Ganopolski, A., Mangini, A., Kubatzki, C., Roth, K., and Kromer, B.: Possible solar origin of the 1,470-year glacial climate cycle demonstrated in a coupled model, Nature, 438, 208–211,, 2005. 

Braun, H., Ditlevsen, P., and Kurths, J.: New measures of multimodality for the detection of a ghost stochastic resonance, Chaos, 19, 1–12,, 2009. 

Brook, E. J., Sowers, T., and Orchardo, J.: Rapid variations in atmospheric methane concentration during the past 110,000 Years, Science, 273, 1087–1091,, 1996. 

Corrick, E. C., Drysdale, R. N., Hellstrom, J. C., Capron, E., Rasmussen, S. O., Zhang, X., Fleitmann, D., Couchoud, I., and Wolff, E.: Synchronous timing of abrupt climate changes during the last glacial period, Science, 369, 963–969,, 2020. 

Crick, L., Burke, A., Hutchison, W., Kohno, M., Moore, K. A., Savarino, J., Doyle, E. A., Mahony, S., Kipfstuhl, S., Rae, J. W. B., Steele, R. C. J., Sparks, R. S. J., and Wolff, E. W.: New insights into the ∼74 ka Toba eruption from sulfur isotopes of polar ice cores, Clim. Past, 17, 2119–2137,, 2021. 

Dansgaard, W., Johnsen, S. J., Clausen, H. B., Dahl-Jensen, D., Gundestrup, N., Hammer, C. U., and Oeschger, H.: North Atlantic climatic oscillations revealed by deep Greenland ice cores, in: Climate processes and climate sensitivity, edited by: Hansen, J. E. and Takahashi, T., American Geophysical Union, Washington D.C., 288–298,, 1984. 

Dansgaard, W., Johnsen, S. J., Clausen, H. B., Dahl-Jensen, D., Gundestrup, N. S., Hammer, C. U., Hvidberg, C. S., Steffensen, J. P., Sveinbjornsdottir, A. E., Jouzel, J., and Bond, G.: Evidence for general instability of past climate from a 250-Kyr ice-core record, Nature, 364, 218–220,, 1993. 

de Bakker, A. T. M., Herbers, T. H. C., Smit, P. B., Tissier, M. F. S., and Ruessink, B. G.: Nonlinear infragravity-wave interactions on a gently sloping laboratory beach, J. Phys. Oceanogr., 45, 589–605,, 2015. 

de Bakker, A. T. M., Tissier, M. F. S., and Ruessink, B. G.: Beach steepness effects on nonlinear infragravity-wave interactions: A numerical study, J. Geophys. Res.-Oceans, 121, 554–570,, 2016. 

Ditlevsen, P. and Crucifix, M.: On the importance of centennial variability for ice ages, Pages Magazine, 3, 152–153,, 2017. 

Ditlevsen, P. D., Svensmark, H., and Johnsen, S.: Contrasting atmospheric and climate dynamics of the last-glacial and Holocene periods, Nature, 379, 810–812,, 1996. 

Ditlevsen, P. D., Kristensen, M. S., and Andersen, K. K.: The recurrence time of Dansgaard-Oeschger events and limits on the possible periodic component, J. Climate, 18, 2594–2603,, 2005. 

Ditlevsen, P. D., Andersen, K. K., and Svensson, A.: The DO-climate events are probably noise induced: statistical investigation of the claimed 1470 years cycle, Clim. Past, 3, 129–134,, 2007. 

Dokken, T. M., Nisancioglu, K. H., Li, C., Battisti, D. S., and Kissel, C.: Dansgaard-Oeschger cycles: Interactions between ocean and sea ice intrinsic to the Nordic seas, Paleoceanography, 28, 491–502,, 2013. 

Elliot, M., Labeyrie, L., and Duplessy, J.-C.: Changes in North Atlantic deep-water formation associated with the Dansgaard–Oeschger temperature oscillations (60–10 ka), Quaternary Sci. Rev., 21, 1153–1165,, 2002. 

Gottwald, G. A.: A model for Dansgaard–Oeschger events and millennial-scale abrupt climate change without external forcing, Clim. Dynam., 56, 227–243,, 2021. 

Griem, L., Voelker, A. H. L., Berben, S. M. P., Dokken, T. M., and Jansen, E.: Insolation and Glacial Meltwater Influence on Sea-Ice and Circulation Variability in the Northeastern Labrador Sea During the Last Glacial Period, Paleoceanogr. Paleocl., 34, 1689–1709,, 2019. 

Hagelberg, T., Pisias, N., and Elgar, S.: Linear and nonlinear couplings between orbital forcing and the marine δ18O record during the late Neogene, Paleoceanography, 6, 729–746,, 1991. 

Hagelberg, T. K., Bond, G., and deMenocal, P.: Milankovitch Band Forcing of Sub-Milankovitch Climate Variability during the Pleistocene, Paleoceanography, 9, 545–558,, 1994. 

Hasselmann, K., Munk, W., and MacDonald, G.: Bispectra of ocean waves, in: Proceedings of the Symposium on Time Series Analysis, edited by: Rosenblatt, M., John Wiley, New York, London, 125–139, 1963. 

Hays, J. D., Imbrie, J., and Shackleton, N. J.: Variations in the Earth's orbit: Pacemaker of the ice ages, Science, 194, 1121–1132,, 1976. 

Herbers, T. H. C. and Burton, M. C.: Nonlinear shoaling of directionally spread waves on a beach, J. Geophys. Res.-Oceans, 102, 21101–21114,, 1997. 

Herbers, T. H. C., Russnogle, N. R., and Elgar, S.: Spectral energy balance of breaking waves within the surf zone, J. Phys. Oceanogr., 30, 2723–2737,<2723:SEBOBW>2.0.CO;2, 2000. 

Hinnov, L. A., Schulz, M., and Yiou, P.: Interhemispheric space-time attributes of the Dansgaard-Oeschger oscillations between 100 and 0 ka, Quaternary Sci. Rev., 21, 1213–1228,, 2002. 

Hodell, D. A., Crowhurst, S. J., Lourens, L., Margari, V., Nicolson, J., Rolfe, J. E., Skinner, L. C., Thomas, N. C., Tzedakis, P. C., Mleneck-Vautravers, M. J., and Wolff, E. W.: A 1.5-million-year record of orbital and millennial climate variability in the North Atlantic, Clim. Past, 19, 607–636,, 2023. 

Huybers, P. and Curry, W.: Links between annual, Milankovitch and continuum temperature variability, Nature, 441, 329–332,, 2006. 

Kelsey, A.: Abrupt climate change and millennial-scale cycles: an astronomical mechanism, Clim. Past Discuss. [preprint],, 2022. 

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,, 2014. 

King, T.: Quantifying nonlinearity and geometry in time series of climate, Quaternary Sci. Rev., 15, 247–266,, 1996. 

Kleppin, H., Jochum, M., Otto-Bliesner, B., Shields, C. A., and Yeager, S.: Stochastic Atmospheric Forcing as a Cause of Greenland Climate Transitions, J. Climate, 28, 7741–7763,, 2015. 

Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A. C. M., and Levrard, B.: A long-term numerical solution for the insolation quantities of the Earth, Astron. Astrophys., 428, 261–285,, 2004. 

Le Treut, H. and Ghil, M.: Orbital forcing, climatic interactions, and glaciation cycles, J. Geophys. Res., 88, 5167–5190,, 1983. 

Liebrand, D. and de Bakker, A. T. M.: Bispectra of climate cycles show how ice ages are fuelled, Clim. Past, 15, 1959–1983,, 2019. 

Lin, J., Svensson, A., Hvidberg, C. S., Lohmann, J., Kristiansen, S., Dahl-Jensen, D., Steffensen, J. P., Rasmussen, S. O., Cook, E., Kjær, H. A., Vinther, B. M., Fischer, H., Stocker, T., Sigl, M., Bigler, M., Severi, M., Traversi, R., and Mulvaney, R.: Magnitude, frequency and climate forcing of global volcanism during the last glacial period as seen in Greenland and Antarctic ice cores (60–9 ka), Clim. Past, 18, 485–506,, 2022. 

Lohmann, G., Butzin, M., Eissner, N., Shi, X., and Stepanek, C.: Abrupt climate and weather changes across time scales, Paleoceanogr. Paleocl., 35, e2019PA003782,, 2020. 

Lohmann, J. and Svensson, A.: Ice core evidence for major volcanic eruptions at the onset of Dansgaard–Oeschger warming events, Clim. Past, 18, 2021–2043,, 2022. 

Loulergue, L., Schilt, A., Spahni, R., Masson-Delmotte, V., Blunier, T., Lemieux, B., Barnola, J.-M., Raynaud, D., Stocker, T. F., and Chappellaz, J.: Orbital and millennial-scale features of atmospheric CH4 over the past 800,000 years, Nature, 453, 383–386,, 2008. 

MacAyeal, D. R.: Binge/purge oscillations of the Laurentide ice sheet as a cause of the North Atlantic's Heinrich events, Paleoceanography, 8, 775–784,, 1993. 

Markle, B. R., Steig, E. J., Buizert, C., Schoenemann, S. W., Bitz, C. M., Fudge, T. J., Pedro, J. B., Ding, Q., Jones, T. R., White, J. W. C., and Sowers, T.: Global atmospheric teleconnections during Dansgaard–Oeschger events, Nat. Geosci., 10, 36–40,, 2017. 

Matyasovszky, I.: Trends, time-varying and nonlinear time series models for NGRIP and VOSTOK paleoclimate data, Theor. Appl. Climatol., 101, 433–443,, 2010. 

Menviel, L. C., Skinner, L. C., Tarasov, L., and Tzedakis, P. C.: An ice-climate oscillatory framework for Dansgaard–Oeschger cycles, Nature Reviews Earth & Environment, 1, 677–693,, 2020. 

Mitsui, T., Lenoir, G., and Crucifix, M.: Is the glacial climate scale invariant?, Dynamics and Statistics of the Climate System, 3, dzy011,, 2019. 

Paillard, D., Labeyrie, L., and Yiou, P.: AnalySeries, Macintosh program performs time-series analysis, EOS Transactions AGU, 77, 379,, 1996. 

Paine, A. R., Wadsworth, F. B., and Baldini, J. U. L.: Supereruption doublet at a climate transition, Communications Earth & Environment, 2, 219,, 2021. 

Pelletier, J. D.: The power spectral density of atmospheric temperature from time scales of 10–2 to 106 yr, Earth Planet. Sc. Lett., 158, 157–164,, 1998. 

Peristykh, A. N. and Damon, P. E.: Persistence of the Gleissberg 88-year solar cycle over the last  12,000 years: Evidence from cosmogenic isotopes, J. Geophys. Res., 108, 1003,, 2003. 

Pestiaux, P., Van der Mersch, I., Berger, A., and Duplessy, J. C.: Paleoclimatic variability at frequencies ranging from 1 cycle per 10 000 years to 1 cycle per 1000 years: evidence for nonlinear behaviour of the climate system, Climatic Change, 12, 9–37,, 1988. 

Petersen, S. V., Schrag, D. P., and Clark, P. U.: A new mechanism for Dansgaard-Oeschger cycles, Paleoceanography, 28, 24–30,, 2013. 

Rasmussen, S. O., Andersen, K. K., Svensson, A. M., Steffensen, J. P., Vinther, B. M., Clausen, H. B., Siggaard-Andersen, M.-L., Johnsen, S. J., Larsen, L. B., Dahl-Jensen, D., Bigler, M., Röthlisberger, R., Fischer, H., Goto-Azuma, K., Hansen, M. E., and Ruth, U.: A new Greenland ice core chronology for the last glacial termination, J. Geophys. Res., 111, 1–16,, 2006. 

Rasmussen, S. O., Bigler, M., Blockley, S. P., Blunier, T., Buchardt, S. L., Clausen, H. B., Cvijanovic, I., Dahl-Jensen, D., Johnsen, S. J., Fischer, H., Gkinis, V., Guillevic, M., Hoek, W. Z., Lowe, J. J., Pedro, J. B., Popp, T., Seierstad, I. K., Steffensen, J. P., Svensson, A. M., Vallelonga, P., Vinther, B. M., Walker, M. J. C., Wheatley, J. J., and Winstrup, M.: A stratigraphic framework for abrupt climatic changes during the Last Glacial period based on three synchronized Greenland ice-core records: refining and extending the INTIMATE event stratigraphy, Quaternary Sci. Rev., 106, 14–28,, 2014. 

Rial, J. A. and Saha, R.: Modeling abrupt climate change as the interaction between sea ice extent and mean ocean temperature under orbital insolation forcing, in: Abrupt Climate Change: Mechanisms, Patterns, and Impacts, Geoph. Monog. Series, 57–74, 193,, 2011. 

Riechers, K., Mitsui, T., Boers, N., and Ghil, M.: Orbital insolation variations, intrinsic climate variability, and Quaternary glaciations, Clim. Past, 18, 863–893,, 2022. 

Rousseau, D.-D., Bagniewski, W., and Ghil, M.: Abrupt climate changes and the astronomical theory: are they related?, Clim. Past, 18, 249–271,, 2022. 

Scafetta, N., Milani, F., Bianchini, A., and Ortolani, S.: On the astronomical origin of the Hallstatt oscillation found in radiocarbon and climate records throughout the Holocene, Earth-Sci. Rev., 162, 24–43,, 2016. 

Schmidt, M. W., Vautravers, M. J., and Spero, H. J.: Rapid subtropical North Atlantic salinity oscillations across Dansgaard–Oeschger cycles, Nature, 443, 561–564,, 2006. 

Schmidt, O. T.: Bispectral mode decomposition of nonlinear flows, Nonlinear Dynam., 102, 2479–2501,, 2020. 

Schulz, M.: On the 1470-year pacing of Dansgaard-Oeschger warm events, Paleoceanography, 17, 1014,, 2002. 

Schulz, M., Berger, W. H., Sarnthein, M., and Grootes, P. M.: Amplitude variations of 1470-year climate oscillations during the last 100,000 years linked to fluctuations of continental ice mass, Geophys. Res. Lett., 26, 3385–3388,, 1999. 

Schulz, M., Paul, A., and Timmermann, A.: Relaxation oscillators in concert: A framework for climate change at millennial timescales during the late Pleistocene, Geophys. Res. Lett., 29, 2193,, 2002. 

Shackleton, N. J., Hall, M. A., and Vincent, E.: Phase relationships between millennial-scale events 64,000-24,000 years ago, Paleoceanography, 15, 565–569,, 2000. 

Stuiver, M. and Braziunas, T. F.: Sun, ocean, climate and atmospheric 14CO2: an evaluation of causal and spectral relationships, The Holocene, 3, 289–305,, 1993. 

Sun, Y., McManus, J. F., Clemens, S. C., Zhang, X., Vogel, H., Hodell, D. A., Guo, F., Wang, T., Liu, X., and An, Z.: Persistent orbital influence on millennial climate variability through the Pleistocene, Nat. Geosci., 14, 812–818,, 2021.  

Svensson, A., Andersen, K. K., Bigler, M., Clausen, H. B., Dahl-Jensen, D., Davies, S. M., Johnsen, S. J., Muscheler, R., Parrenin, F., Rasmussen, S. O., Röthlisberger, R., Seierstad, I., Steffensen, J. P., and Vinther, B. M.: A 60 000 year Greenland stratigraphic ice core chronology, Clim. Past, 4, 47–57,, 2008. 

University of Copenhagen, Centre for Ice and Climate, Niels Bohr Institute:,,, last access: 15 July 2023. 

van Hoesel, A., Hoek, W. Z., Pennock, G. M., and Drury, M. R.: The Younger Dryas impact hypothesis: a critical review, Quaternary Sci. Rev., 83, 95–114,, 2014. 

van Kreveld, S., Sarnthein, M., Erlenkeuser, H., Grootes, P., Jung, S., Nadeau, M. J., Pflaumann, U., and Voelker, A.: Potential links between surging ice sheets, circulation changes, and the Dansgaard-Oeschger Cycles in the Irminger Sea, 60–18 Kyr, Paleoceanography, 15, 425–442,, 2000. 

Vettoretti, G. and Peltier, W. R.: Fast Physics and Slow Physics in the Nonlinear Dansgaard–Oeschger Relaxation Oscillation, J. Climate., 31, 3423–3449,, 2018. 

Vettoretti, G., Ditlevsen, P., Jochum, M., and Rasmussen, S. O.: Atmospheric CO2 control of spontaneous millennial-scale ice age climate oscillations, Nat. Geosci., 15, 300–306,, 2022. 

Vinther, B. M., Clausen, H. B., Johnsen, S. J., Rasmussen, S. O., Andersen, K. K., Buchardt, S. L., Dahl-Jensen, D., Seierstad, I. K., Siggaard-Andersen, M.-L., Steffensen, J. P., Svensson, A., Olsen, J., and Heinemeier, J.: A synchronized dating of three Greenland ice cores throughout the Holocene, J. Geophys. Res., 111, D13102,, 2006. 

Wagner, G., Beer, J., Masarik, J., Muscheler, R., Kubik, P. W., Mende, W., Laj, C., Raisbeck, G. M., and Yiou, F.: Presence of the solar de Vries cycle (205 years) during the last ice age, Geophys. Res. Lett., 28, 303–306,, 2001. 

Wara, M. W., Ravelo, A. C., and Revenaugh, J. S.: The pacemaker always rings twice, Paleoceanography, 15, 616–624,, 2000. 

Wolff, E. W., Chappellaz, J., Blunier, T., Rasmussen, S. O., and Svensson, A.: Millennial-scale variability during the last glacial: The ice core record, Quaternary Sci. Rev., 29, 2828–2838,, 2010. 

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. 

Short summary
Climate cycles with millennial periodicities are enigmatic because no Earth external climate forcing exists that operates on millennial timescales. Using a statistical analysis of a famous Greenlandic air temperature record, we show that two disparate energy sources (one astronomical and one centennial) fuel millennial climate variability. We speculate that two distinct Earth internal cryospheric/climatic/oceanic processes are responsible for the transfer of energy to millennial climate cycles.