Abrupt climate changes and the astronomical theory: are they related?

​​​​​​​Abrupt climate changes are defined as sudden climate changes that took place over tens to hundreds of years or recurred at millennial timescales; they are thought to involve processes that are internal to the climate system. By contrast, astronomically forced climate changes involve processes that are external to the climate system and whose multi-millennial quasi-periodic variations are well known from astronomical theory. In this paper, we re-examine the main climate variations determined from the U1308 North Atlantic marine record, which yields a detailed calving history of the Northern Hemisphere ice sheets over the past 3.2 Myr. The magnitude and periodicity of the ice-rafted debris (IRD) events observed in the U1308 record allow one to determine the timing of several abrupt climate changes, the larger ones corresponding to the massive iceberg discharges labeled Heinrich events (HEs). In parallel, abrupt warmings, called Dansgaard–Oeschger (DO) events, have been identified in the Greenland records of the last glaciation cycle. Combining the HE and DO observations, we study a complex mechanism giving rise to the observed millennial-scale variability that subsumes the abrupt climate changes of last 0.9 Myr. This process is characterized by the presence of Bond cycles, which group DO events and the associated Greenland stadials into a trend of increased cooling, with IRD events embedded into every stadial, the latest of these being an HE. These Bond cycles may have occurred during the last 0.9 Ma when Northern Hemisphere ice sheets reached their maximum extent and volume, thus becoming a major player in this time interval's climate dynamics. Since the waxing and waning of ice sheets during the Quaternary period are orbitally paced, we conclude that the abrupt climate changes observed during the Middle Pleistocene and Upper Pleistocene are therewith indirectly linked to the astronomical theory of climate.

1 Introduction

Well-dated geological data indicate that the Earth experienced orbitally paced climate changes from at least the late Precambrian, i.e., 1.4 billion years ago during the Proterozoic Eon (Benn et al., 2015; Zhang et al., 2015; Hoffman et al., 2017; Meyers and Malinverno, 2018), and all along the Phanerozoic (Lisiecki and Raymo, 2005; Liebrand et al., 2011; Miller et al., 2011; Kent et al., 2017, 2018; Olsen et al., 2019; Drury et al., 2021; Westerhold et al., 2020). These changes reflect the variations in the Earth's axis of rotation – namely in its precession and tilt – and in the geometry of the Earth's orbit around the sun, i.e., in its eccentricity, driven by gravitational interactions within the solar system (Berger, 1977; Laskar et al., 2011; Hinnov, 2013, 2018). These variations affect the distribution of insolation at the top of the atmosphere, forcing latitudinal and seasonal climate changes with periodicities of tens or hundreds of thousands of years. Although it is presently well acknowledged in the climate community, this mechanism underwent a long journey to reach such general acceptance (Imbrie and Imbrie, 1979).

In the 1840s, while visiting Scotland, Agassiz (1842) attributed erratic boulders noticed in the local landscape to former ice age glaciers in the area, recalling his former observations in Switzerland (Agassiz, 1838). Meanwhile, Adhémar (1842) proposed that glaciations, like those inferred by Agassiz from field observations, occurred every 22 000 years due to the Earth's evolving on an elliptic orbit with its rotation axis tilted with respect to the orbital plane. Adhémar also concluded that, as the Southern Hemisphere receives less solar radiation per year than the Northern Hemisphere (Adhemar was relying on the number of the nights at the boreal pole versus the number of nights at the austral pole), this contributed to keeping temperatures cold enough to allow ice sheets to build up.

It was 5 decades later that Croll (1890) presented his theory of the ice ages being driven by the changing distance between the Earth and the Sun as measured on 21 December, which is due to the eccentricity of the Earth's orbit and the precession of the equinoxes. He used formulae for orbital variations developed by Le Verrier (1858) and also stated that, above a particular threshold, Northern Hemisphere winters would trigger an ice age, while below another threshold an ice age would develop in the Southern Hemisphere. According to Croll's theory, it is the slowly evolving eccentricity of the Earth's orbit that is the key driver, with glaciations occurring only when eccentricity is high. He concluded, therefore, that eccentricity can impact the annual amount of heat received from the sun and also leads to differences in seasonal temperatures.

Parallel to these physical calculations, observations of frontal moraines in the Alpine valleys led to the identification of numerous glacial events in the past. Considering the geometry and position of the terrace systems resulting from these moraines in Alpine valleys, Penck and Brückner (1909) identified four main glaciations, determined from ice advances in the Alpine foreland: Günz, Mindel, Riss, and Würm. This sequence became the paleoclimate framework for many decades, until the development of marine and ice coring programs.

A total of 30 years after Croll, Milankovitch (1920, 1941), inspired by his exchanges with Wladimir Peter Köppen and Alfred Wegener​​​​​​​, provided a crucial step in advancing the theory that major past climate changes had an astronomical origin resulting from the interaction between the eccentricity of the Earth's orbit around the sun, the precession of the equinoxes, and the tilt of the Earth's rotation axis. This theory, often given his name today, implies a fairly regular, multi-thousand-year variability. Based on this crucial idea and on using the recent orbital calculations of Pilgrim (1904), among others, Milankovitch was able to estimate temperatures and insolation at various latitudes at the top of the atmosphere and follow their variations through time. Among the results of these estimates, he found that the summer insolation at 65^∘ N is best correlated with glacial–interglacial transitions, as determined by the study by Penck and Brückner (1909) of Alpine glaciations.

Milankovitch's astronomical theory of climate was severely criticized by contemporary physicists and rejected by most Quaternary geologists, as isotopic C¹⁴ dating threw the Günz–Mindel–Riss–Würm classification out the window (e.g., Flint, 1971, and references therein). Interest in this theory was restored, however, as the pace of Quaternary climate changes – defined now by δ¹⁸O records from benthic foraminifera (Emiliani, 1955) – was connected on a much sounder basis with orbital variations in the seminal paper of Hays et al. (1976). At the same time, Berger (1977, 1978), using the much more accurate orbital calculations of Chapront et al. (1975), more reliably linked the insolation variations to the waxing and waning of the huge Northern Hemisphere ice sheets over North America, Greenland, Iceland, and Europe, including the British Isles and Fennoscandia. The Quaternary period, however, also shows an intriguing transition between the 40 kyr cycles that is dominant during the Early (or Lower) Pleistocene, i.e., from 2.6 up to 1.25 Ma, and the dominant 100 kyr cycles of the Middle Pleistocene and Late (or Upper) Pleistocene, i.e., the last 800 kyr, which show longer glacials that imply much larger continental ice sheets. The transition interval between 1.25 and 0.8 Ma, named the mid-Pleistocene transition (MPT; Pisias and Moore, 1981; Ruddiman et al., 1989), experienced variability that does not appear to be directly linked to insolation changes at high latitudes. In fact, the processes involved are not fully understood and are still a matter of debate (Ghil and Childress, 1987, Sect. 12.7; Ghil, 1994; Saltzman, 2002; Clark et al., 2021).

Although the broad astronomic framework for past climate changes seems to be widely accepted, high-resolution investigations over the past decades in ice, marine, and terrestrial records revealed much shorter periodicities than the orbital ones, as well as abrupt changes that do not match the transitions attributed to variations in eccentricity, obliquity, or precession of the equinoxes. This millennial variability, identified initially during the last climate cycle (Dansgaard et al., 1982), has been described in older records, indicating that it prevailed during at least the past 0.5 Ma, as recorded in marine records (McManus et al., 1999). The increased resolution of these records has shown that millennial climate variability can be observed in all types of records, including in ice (Johnsen et al., 2001; Buizert et al., 2015a), marine (Broecker, 2000; Shackleton et al., 2000), and terrestrial (Sanchez-Goni et al., 2000, 2002; Cheng et al., 2016) records associated with abrupt warmings or with strong (and massive) iceberg discharges into the North Atlantic Ocean (Ruddiman, 1977; Bond et al., 1992, 1993; Andrews and Voelker, 2018) that correspond to particular cooling events. Both hemispheres have been impacted by this millennial variability through a mechanism described as the bipolar seesaw (Stocker and Johnsen, 2003).

In this paper, we show that abrupt climate changes are still affected, albeit indirectly, by changes in insolation and hence in overall ice sheet volume. In this sense, the Milankovitch framework is shown herein to be quite relevant to the study of abrupt and large climate changes. The paper is organized as follows. In Sect. 2, we describe the records and the methods used to identify the major transitions. In Sect. 3, major transitions for the last 3.2 Myr of Northern Hemisphere climate are described. In Sect. 4, we concentrate on the millennial-scale variability revealed by these records. In Sect. 5, we establish a connection between Dansgaard–Oeschger (DO) events and amended Bond cycles, and outline how global ice sheet volume affects the latter. Concluding remarks follow in Sect. 6.

2 Material and methods

2.1 The records

Over the past 66 Myr, corresponding to the Cenozoic Era, Earth's climate has experienced four main states, from “Warmhouse” and “Hothouse” between 66 and 34 Ma to “Coolhouse” and “Icehouse” from 34 Ma until the present time; see Fig. 1a. Although the first two states alternated in a warm–hot–warm sequence, the last two succeeded each other, thus generating the classical climate trend towards the recent ice age conditions (Zachos et al., 2001; Westerhold et al., 2020; Scotese et al., 2021). The last 3.3 Myr have been defined as an Icehouse climate state, with the appearance of the Northern Hemisphere ice sheets and their variations through time (Westerhold et al., 2020). This Icehouse state is characterized by a change of the interplay between benthic δ¹³C and δ¹⁸O, which corresponds to a new relationship between the carbon cycle and climate. Indeed, the negative excursions in δ¹³C were associated with negative ones in δ¹⁸O during most of the Cenozoic since 66 Ma. A shift occurred, however, in the Plio-Pleistocene at about 5 Ma, with negative δ¹³C excursions being associated with positive δ¹⁸O excursions (Turner, 2014). Such a shift appears to be related to a dichotomy in the response of the marine and terrestrial carbon reservoirs to orbital forcings.

Figure 1 Earth climate history represented by two marine records and an ice core one. (a) Record of benthic δ¹⁸O of the CENOGRID stack for the last 66 Myr; the hothouse, warmhouse, coolhouse, and icehouse intervals, following Westerhold et al. (2020), are indicated along the abscissa. (b) Benthic δ¹⁸O record of the U1308 marine core (Hodell and Channell, 2016a) for the last 3.2 Myr; the 40 kyr, 100 kyr, and mid-Pleistocene transition (MPT) intervals are indicated along the top of the diagram. (c) Record of 122 ka b2k δ¹⁸O of the NGRIP ice core (Rasmussen et al., 2014) for the last 122 kyr b2k. Here “b2k” means before 2000 CE.

The past 3.2 Myr of Northern Hemisphere climate are particularly well described in North Atlantic core U1308 at 49.87^∘ N, 24.24^∘ W (Hodell and Channell, 2016a). This core is located in the ice-rafted debris (IRD) of continental detrital material eroded by the ice sheets (Ruddiman, 1977), and it yields a more complete record than U1313 (Naafs et al., 2013); see Fig. 1b. The variations in the benthic δ¹⁸O mostly indicate varying periodicities through time that correspond to periodicities in the orbital parameters of Earth's climate (Hodell and Channell, 2016a; see Fig. S1 in the Supplement), as confirmed by Lisiecki and Raymo (2005) from the stack oxygen isotope record that they produced using 57 marine records from the world's oceans.

The behavior of the U1308 proxy records on the timescale of many tens and hundreds of thousands of years allows us to track the numerous glacial–interglacial cycles of the past 2.75 Myr with the help of variations in Earth's orbital parameters (Fig. S1). Besides these relatively slow variations, evidence of millennial-scale variability can be observed since the appearance of IRD in the North Atlantic at about 1.5 Ma. This much faster variability is superposed upon the classical orbital periodicities; see Fig. 3.

Observations of such abrupt variations have been reported in some detail for the last glacial period by the North Greenland Ice Core Project at 75.1^∘ N, 42.32^∘ W (Barbante et al., 2006). The δ¹⁸O record best describes the more or less regular recurrence of cold and warm events that we are analyzing; see Fig. 1c.

2.2 The methods

The increase in resolution of the Greenland ice cores (Fischer et al., 2015; Schupbach et al., 2018; Svensson et al., 2020) and of several speleothems (Cheng et al., 2016) allowed much more detailed analyses not only of the δ¹⁸O but also of other proxies, such as dust content, leading to the identification of sub-events, as compiled by Rasmussen et al. (2014). It is these high-resolution analyses that led to defining the Greenland interstadials (GIs), followed by the associated Greenland stadials (GSs). To gain further insight into the climate story the records tell us, we performed a quantitative, objective analysis of these time series of proxy variables based on two methods: the Kolmogorov–Smirnov (KS) augmented test of Bagniewski et al. (2021) and the recurrence plots (RPs) of Marwan et al. (2007, 2013).

Bagniewski et al. (2021) proposed a method based on an augmented nonparametric Kolmogorov–Smirnov (KS) test to efficiently and robustly detect abrupt transitions in paleo-records. The classical KS test consists of quantifying and comparing the empirical distribution functions of two samples from a time series both before and after a potential jump. In the present case, the test is augmented by varying window sizes and by evaluating the rate of change and the trend in maxima and minima of the time series. This method allows one to reliably establish the main transitions in a record, such as the marine isotope stages in the U 1308 benthic δ¹⁸O or the GS–GI boundaries in the North Greenland Ice Core Project (Barbante et al., 2006) δ¹⁸O.

The recurrence plots (RPs) were introduced by Eckmann et al. (1987) into the study of dynamical systems and popularized in the climate sciences by Marwan et al. (2007, 2013). The purpose of RPs is to identify recurring patterns in a time series in general and in a paleoclimate time series in particular.

The RP for a time series $x_i:i=\mathrm{1}$, …, N is constructed as a square matrix in a Cartesian plane with the abscissa and ordinate both corresponding to a time-like axis, with one copy x_i of the series on the abscissa and another copy x_j on the ordinate. A dot is entered into a position (i,j) of the matrix when x_j is sufficiently close to x_i. For more details – such as how “sufficiently close” is determined – we refer to Eckmann et al. (1987) and to Marwan et al. (2013). Clearly, all the points on the diagonal i=j have dots and the matrix is in general rather symmetric, although one does not always define closeness symmetrically; to wit, x_j may be “closer to” x_i than x_i is to x_j (Eckmann et al., 1987). An important advantage of the RP method is that it does apply to dynamical systems that are not autonomous, i.e., that may be subject to time-dependent forcing. The latter is certainly the case for the climate system on timescales of 10–100 kyr and longer, which are affected strongly by orbital forcing.

Eckmann et al. (1987) distinguished between large-scale “typology” and small-scale “texture” in the interpretation of square matrix of dots that is the visual result of RP. Thus, if all the characteristic times of an autonomous dynamical system are short compared to the length of the time series, the RP's typology will be homogeneous and thus not very interesting. In the presence of an imposed drift, a more interesting typology will appear. The most interesting typology in RP applications so far is associated with recurrent patterns that are not quite (but are close to being) periodic. Hence, such patterns are not easily detectable by purely spectral approaches to time series analysis. Marwan et al. (2013) discuss how to render the purely visual RP typologies studied up to that point in a more objectively quantifiable way via recurrence quantification analysis and bootstrapping (Efron, 1981; Efron and Tibshirani, 1986).

Table 1 Statistics of the recurrence rate analysis for U1308 benthic δ¹⁸O, U1308 bulk carbonate δ¹⁸O (Hodell and Channell, 2016a), and NGRIP δ¹⁸O (Rasmussen et al., 2014). Dates where the recurrence rate (RR) prominence is higher than the RR standard deviation are considered major transitions and are labeled with a red cross in Figs. 2c, 3c, and 4c. Periods with an RR prominence that is slightly lower than the RR standard deviation are taken into account and labeled with a pink cross in the same figures. RR prominence significance is represented in this table as follows: ^a shows prominence significance above the standard deviation (SD), and ^b shows prominence significance slightly lower than the standard deviation.

The selection of the major transitions in the RPs was achieved by performing a recurrence rate (RR) analysis. The values obtained for every record correspond to the mean for different window lengths ranging from 1 to 15 kyr. The selection of the transitions of interest thus relies on the definition of a threshold that we choose to be the standard deviation of RR prominence, which is 0.089 for the U1308 benthic δ¹⁸O, 0.127 for the U1308 bulk carbonate δ¹⁸O, and 0.173 for NGRIP δ¹⁸O. The results were plotted along with the original recurrence plot (Table 1). Please see further details in Bagniewski et al. (2021).

3 The past 3.2 Myr history of the Northern Hemisphere climate

Analyzing the U1308 benthic δ¹⁸O record allows one to identify the various transitions between the marine isotope stages that covered the last 3.2 Ma and that were driven by variations in the orbital parameters (Fig. S1). However, the RPs reveal several key features that appear in the U1308 records during this time interval that are much more plausibly related to processes other than the orbital forcing. Hoddell and Channell (2016a) identified four major steps in their climate record. These four steps are linked to thresholds in the benthic and bulk carbonate δ¹⁸O variations, and they occur at 2.75, 1.5, 0.9, and 0.65 Ma.

The first date is interpreted as corresponding to the earliest occurrence of IRD in the North Atlantic. This occurrence characterizes the presence of Northern Hemisphere coastal glaciers large enough to calve icebergs into the ocean, and the melting of these icebergs is likely to have impacted the oceanic circulation. Naafs et al. (2013) reported the occurrence of minor IRD events, attributed mainly to Greenland and Fennoscandian glaciers. These events indicate a greater size of ice sheets over these regions than during the later Quaternary, when the North American ice sheets were considerably larger.

The second date corresponds to an increase in amplitude in ice volume variations between glacial maxima and interglacial optima. This second step shows the permanent occurrence of ice-rafted events during glacial intervals in the record. The latter behavior signals a stronger relationship between climate variations and Northern Hemisphere ice sheets.

The third date, close to the MIS22–24 δ¹⁸O optima, shows increased continental ice volume in the Northern Hemisphere (Batchelor et al., 2019) but also more stability in the East Antarctic ice sheet in Southern Hemisphere (Jakob et al., 2020). In parallel, evidence of a major glacial pulse recorded in Italy's Po Plain, as well as in ¹⁰Be-dated boulders in Switzerland, is interpreted as marking the onset of the first major glaciation in the Alps (Muttoni et al., 2003; Knudsen et al., 2020). At about the same time, the synthetic Greenland δ¹⁸O reconstruction of Barker et al. (2011), which starts at 800 ka, indicates the occurrence of millennial variability expressed by DO-like events.

The last date at 0.65 Ma marks the end of the transition from the Lower Pleistocene and Middle Pleistocene interval – characterized by 41 kyr dominated cycles and smaller 23 kyr ones – to the Upper Pleistocene, with its 100 kyr dominated cycles; see Fig. 1b. The sawtooth pattern of the interglacial–glacial cycles (Broecker and van Donk, 1970), which first becomes noticeable at 0.9 Ma, is well established during this final interval, in contradistinction with the previous, more smoothly shaped pattern that appears to follow the obliquity variations. The global ice volume is maximal, exceeding the values observed earlier in the record, due in large part to the larger contribution of the Northern American ice sheets. The latter now have a bigger impact on Northern Hemisphere climate than the Eurasian ice sheets (Batchelor et al., 2019). The IRD event intensity and frequency of occurrence increased (McManus et al., 1999) as well, leading to the major iceberg discharges into the North Atlantic named Heinrich events (HEs); see Heinrich (1988), Bond et al. (1992, 1993), and Obrochta et al. (2014). The interval of 1–0.4 Ma is also the interval during which Northern Hemisphere ice sheets reached a southernmost extent during MIS16 and MIS12, similar to the one reached during MIS6 (Batchelor et al., 2019).

Table 2 Comparison of the main steps detected by Hodell and Channell (2016a) from the U1308 marine record with those deduced from the recurrence plot of the benthic δ¹⁸O and bulk carbonate δ¹⁸O data of the same record. RR prominence significance is defined the same way as in Table 1.

Figure 2 Recurrence analysis of the δ¹⁸O record in U1308 Cibicidoides sp. (Hodell and Channell, 2016a). (a) Time series of the Cibicidoides sp. δ¹⁸O (blue curve, top) and of the bulk carbonate δ¹⁸O (magenta curve, bottom). (b) Recurrence plot (RP) where the proximity threshold is Δ(δ¹⁸O)=0.2 ‰. (c) Recurrence rate (RR) with significant values labeled with a red cross; RR close to the threshold is marked in pink. Vertical bars represent the major transitions (5 + 1) determined by the analysis. The RP website is http://www.recurrence-plot.tk/ (last access: 8 February​​​​​​​ 2022).​​​​​​​

The benthic δ¹⁸O record of the U1308 marine-sediment core is interpreted in terms of global ice volume and deep-ocean temperatures (Chappell and Shackleton, 1986; Shackleton, 2000; Elderfield et al., 2012). Its recurrence analysis shows a drift topology (Marwan et al., 2007, 2013) that characterizes a monotonic trend in time, associated with nonstationary systems with slowly varying parameters. Moreover, the RP exhibits a characteristic texture given by the pattern of vertical and horizontal lines that mark recurrences. These lines sometimes form recurrence clusters that correspond to specific periodic patterns. We thus identify five steps in the δ¹⁸O variability; see Fig. 2a and c and Tables 1 and 2. Two are roughly similar to those determined by Hoddell and Channell (2016a), i.e., at 1.5 and 0.65 Ma, and three differ, i.e., at 2.55, 1.25, and 0.35 Ma. Interestingly, the interval 1.25 to 0.65 Ma corresponds roughly to the previously mentioned MPT, during which a shift from climate cycles dominated by a 40 kyr periodicity to 100 kyr dominated ones occurred (Shackleton and Opdyke, 1977; Pisias and Moore, 1981; Ruddiman et al., 1989; Clark and Pollard, 1998; Clark et al., 2006). A sixth transition at 2.95 Ma is close to the selected threshold, as described at the end of Sect. 2.2.

Figure 3 Recurrence analysis of U1308 bulk carbonate δ¹⁸O (Hodell and Channell, 2016a). Panel (a) is the same as Fig. 2a.(b) RP of the bulk carbonate δ¹⁸O, with the same proximity threshold as in Fig. 2a. Panel (c) is the same as in Fig. 2c, except that the vertical bars now represent the major thresholds (4 + 3) determined by the RR analysis in panel (c). The RP website is http://www.recurrence-plot.tk/.

The δ¹⁸O bulk carbonate record in the U1308 core is in turn interpreted as characterizing IRD released into the North Atlantic Ocean, with the most negative δ¹⁸O values representing the largest iceberg calvings (Hodell and Channell, 2016a). The recurrence analysis of this record also displays the evolutionary trend represented by a drift topology, and it yields one of the two Hodell and Channell (2016a) steps at 1.5 Ma, while the transition at 0.65 Ma is slightly below our threshold (Fig. 3a, c, Tables 1 and 2). Our analysis further identifies the steps at 1.25, 1.65 and 2.75 Ma, with the 1.25 step also noticed in the δ¹⁸O. Two other transitions show a recurrence prominence close to our threshold at 1.95 and 2.45 Ma. Table 2, comparing the Hodell and Channell (2016a) steps with the thresholds detected by the RP, shows that the former are mainly related to the IRD history of the past 3.3 Ma.

The 1.25 Ma date is particularly significant, since it is followed by an increase in the amplitude of glacial–interglacial fluctuations (Fig. 2a). The interval 2.8 to 1.2 Ma shows glacial–interglacial sea level variations of about 25–70 m below the present-day value (van de Wal et al., 2011). The reconstructed CO₂ concentrations varied between 270 and 280 ppmv during interglacials and between 210 and 240 ppmv during glacials, with a decreasing trend of about 23 ppmv over this 1.4 Myr long interval (van de Wal et al., 2011). After 1.25 Ma, the sea level changes increased to about 70–120 m below the present-day value, while the reconstructed CO₂ concentrations varied between 250 and 320 ppmv during interglacials and between 170 and 210 ppmv during glacials (Berends et al., 2021). Similar variations were determined by Seki et al. (2010), although pCO₂ changes that occurred before the time reached by ice core records are associated with high uncertainties in both dating and values. The “Milankovitch glacials”, which correspond to the odd marine isotope stages determined in the U1308 core and in many others, have maxima that are characterized by low eccentricity and obliquity, and a boreal summer that coincides with aphelion and leads and, therefore, to minimum values of summer insolation. The increase in IRD variability and magnitude since 1.5 Ma (Hodell and Channell, 2016a), however, shows that distinct, faster processes have to be considered due to slow changes in Earth's orbital parameters; see again Figs. 2 and 3.

4 Millennial-scale variability

The behavior of the U1308 proxy records on the timescale of many tens and hundreds of thousands of years was described briefly in the previous section, and it allowed us to track the numerous glacial–interglacial cycles of the past 2.75 Myr with the help of variations in Earth's orbital parameters; see Fig. S1 in the Supplement. Besides these relatively slow variations, evidence of millennial-scale variability can be observed since the appearance of IRD in the North Atlantic at about 1.5 Ma; this much faster variability is superposed upon the classical orbital periodicities (Fig. 3).

Detailed observations of such abrupt variations have been reported for the last glacial period, with the more or less regular recurrence of cold and warm events; see Fig. 1c. The former are represented by IRD events, some of which are significantly stronger and represent the previously mentioned HEs that correspond to massive discharges of icebergs into the North Atlantic (Heinrich, 1988; Bond et al., 1992; McManus et al., 1994; Hemming, 2004).

Figure 4 Recurrence analysis of NGRIP δ¹⁸O. (a) NGRIP δ¹⁸O variations over the last 122 kyr b2k (Rasmussen et al., 2014); selected canonical Dansgaard-Oeschger (DO) events are indicated on the abscissa by the numbers assigned by Dansgaard et al. (1993). (b) RP of the time series in panel (a) above; the proximity threshold here is Δ(δ¹⁸O)=1.3 ‰. (c) Recurrence rate analysis showing the selected thresholds (7 + 7) indicated by the vertical bars. The RP web site is http://www.recurrence-plot.tk/.

Abrupt warmings happening over as little as a few decades each have been inferred from Greenland ice cores (Dansgaard et al., 1993; Clark et al., 1999) and labeled DOs or Greenland interstadials (GIs: Rasmussen et al., 2014); see Fig. 4a. These warm events are followed by a return to glacial conditions, called Greenland stadials (GSs). This return generally happens in two steps, thus forming DO cycles of variable duration that do not exceed a millennial timescale (Broecker, 1994; Boers et al., 2018; Boers, 2018, and references therein). Broecker (1994) does not provide any precise DO timescale but indicates that “climate cycles [average] a few thousands of years in duration”. Boers et al. (2018) do not yield any exact value for the DO timescales either, but Fig. S1 in their appendix does indicate the durations of the stadials and interstadials in the NGRIP record varying between 340 and 8365 years and between 190 and 16 440 years respectively.

Table 3 Thresholds identified in the recurrence plot of NGRIP δ¹⁸O record and their correspondence in the marine isotope stratigraphy of the last climate cycle stratigraphy (Bassinot et al., 1994; Lisiecki and Raymo, 2005). RR prominence significance as in Table 1.

Numerous DO timescales have been published by Rasmussen et al. (2014), Wolff et al. (2010), Rousseau et al. (2017a, b), and most recently Capron et al. (2021). However, no DO cycle timescale has been published as of yet. Table 3 reports the duration of the DO cycles as expressed by the interval between the start of a GI and the end of a GS as estimated by Rasmussen et al. (2014). It shows that the average duration of a DO cycle is of 4045 ± 3179 years. As the periodicity of the events is at millennial and sub-millennial scale (Broecker, 1994; Clark et al., 1999; Ganopolski and Rahmstorf, 2001; Rahmstorf, 2002; Schulz, 2002; Menviel et al., 2014; Lohmann and Ditlevsen, 2018, 2019) it corresponds to processes that cannot be related to any orbital forcing (Lohmann et al., 2020) but rather to factors that are intrinsic to the Earth system.

DO events were first observed in the various ice cores retrieved from the Greenland ice sheet (Dansgaard et al., 1969; Johnsen et al., 1972). Some of these DOs were correlated to European warm interstadials, which had been described from pollen records (Woillard, 1978; Behre, 1989; Zagwijn, 1989). The existence and dating of the DOs was initially questioned as they had not been observed in marine cores in the 1970s and 1980s (Broecker et al., 1988; Broecker and Denton, 1989). However, Dansgaard et al. (1993) clearly identified 23 rapid warming events during the last climate cycle from the Greenland GRIP ice core.

These 23 DO events were later confirmed in other Greenland ice cores (Johnsen et al., 2001), following the initial correlation of rapid changes between Camp Century and Dye3 ice cores (Dansgaard et al., 1982), and they are considered the “canonical” DOs. They were assigned numbers that increase sequentially downcore, with no. 1 allocated to the Bølling pollen oscillation (Dansgaard et al., 1993); see Fig. 4a. Moreover, the 23 canonical DOs were described later using various types of records in marine sediments (Bond et al., 1992; Henry et al., 2016) and in terrestrial sediments (Allen et al., 1999; Sanchez-Goni et al., 2000, 2002; Müller et al., 2003; Fletcher et al., 2010; Rousseau et al., 2017a, b, 2021), including speleothems (Wang et al., 2001; Genty et al., 2003; Fleitmann et al., 2009; Boch et al., 2011). The transitions identified by a KS test in the NGRIP δ¹⁸O include all the canonical events described by Dansgaard et al. (1993) and identified in Rasmussen et al. (2014), although a few other DOs or sub-events are not detected by it even after changing the test window size; see Fig. S2.

Using a global climate indicator like methane (CH₄) Blunier and Brook (2001), EPICA Community Members (Barbante et al., 2006) and the WAIS Consortium (Buizert et al., 2015a, b) demonstrated that the millennial-scale variations observed during the last 130 kyr in Greenland are observed in Antarctica as well and that they are thus a global phenomenon. Hinnov et al. (2002) also carried out an investigation on the methane-linked GISP2 and Byrd ice cores, including a statistically constrained spectral coherency analysis demonstrating the global reach of the DO cycles. The δ¹⁸O variations in the two hemispheres are in opposite phases though, with the Southern Hemisphere warmings occurring prior to the Northern Hemisphere ones.

Several hypotheses have been proposed to determine whether the climatic signal propagated between the two hemispheres is in a southward or northward direction. Variations in the Atlantic meridional overturning circulation (AMOC) play a key role in this teleconnection, as AMOC slowdown during GSs corresponds to reduced oceanic northward heat transport (Sarnthein et al., 2001; Ganopolski and Rahmstorf, 2001; McManus et al., 2004). This reduction does not contradict the southward propagation direction of the climate signal, as documented by detailed high-resolution marine-sediment and ice core studies (Buizert et al., 2015b; Henry et al., 2016), with the Greenland climate leading Antarctica by approximately 200 years. On the other hand, Knorr and Lohmann (2003) suggest a true south to north direction of climate change propagation. At this point, one can only state that the latter observational studies do not contradict the earlier modeling studies but also do not directly confirm them either.

Many DO models – e.g., Buizert and Schmittner (2015), Dokken et al. (2013), Ganopolski and Rahmstorf (2001), Lohman and Ditlevsen (2018, 2019), Peltier and Vettoretti (2014), Shaffer et al. (2004), Klockmann et al. (2018), Menviel et al. (2014, 2021), and Timmermann et al. (2003) – have not specifically addressed the issue of the interhemispheric signal's direction. To address this issue, Boers et al. (2018) recently developed a simple model to reconstruct the millennial variability in δ¹⁸O of the past 60 kyr b2k, as observed in the high-resolution ice cores from NGRIP in Greenland and the West Antarctic ice sheet (WAIS) in Antarctica. This simple model, based on the bipolar seesaw mechanism (Stocker and Johnsen, 2003), combines the interactions between ice shelves and sea ice extents, subsurface water temperatures in the North Atlantic Ocean, atmospheric temperature in Greenland, AMOC strength, and δ¹⁸O in both Greenland and Antarctica. The interplay of the feedbacks involved allowed the authors to reproduce the millennial-scale variability observed in both hemispheres despite the lack of any time-dependent forcing that would involve orbital parameters.

5 DO events and bond cycles

The results outlined in the previous section have led us to concentrate on the canonical DOs for which the temperature reconstruction has been proposed based on ¹⁵N measurements from the Greenland ice (Guillevic et al., 2014). These estimates indicate that the 21 warming events during the last climate cycle, between 12 and 87 ka b2k, had a range of 10–12 ^∘C on average (Kindler et al., 2014), with each transition lasting between 50 and 100 years on average (Wolff et al., 2010; Rousseau et al., 2017a, b, 2021), as found at least near the top of the Greenland ice sheet at the coring site. Not only does the change in temperature over such a short time imply a drastic reorganization of the atmospheric and associated marine circulations (Boers et al., 2018), but the timing does not correspond to any periodicity of the orbital parameters, either in the average duration of the events or in the average interval between two such events; see Fig. S2.

Figure 5 Recurrence analysis of NGRIP δ¹⁸O. (a) NGRIP δ¹⁸O variations over the last 122 kyr b2k (Rasmussen et al., 2014); same as in Fig. 4a. Panel (b) is the same as Fig. 4b. (c) Variation of the sea level over the last 122 kyr BP as reconstructed from benthic δ¹⁸O foraminifera by Waelbroeck et al. (2002). The bold red line is the global mean sea level (m) below present, and the light blue lines are the minimum and maximum global mean sea level values (m) below present. The horizontal green line segments indicate some key sea levels. The RP web site is http://www.recurrence-plot.tk/.

The recurrence analysis performed on the high-resolution NGRIP δ¹⁸O record and illustrated in Figs. 4b and 5b suggests looking for other mechanisms that might cause the abrupt changes mentioned previously. As for the δ¹⁸O record analysis in core U1308 displayed in Figs. 2 and 3, the NGRIP study shows a drift topology with a non-uniform pattern. The RP shows changes in the system's regime of behavior, as identified by transitions detected by analyzing the RR values. The RR analysis of the NGRIP δ¹⁸O values identified seven major transitions. Seven more indicate an RR close to the selected threshold; see Tables 1 and 3. Numerous thresholds identified correspond to key dates of the last climate cycle stratigraphy (Bassinot et al., 1994; McManus et al., 1994; Kukla et al., 1997; Lisiecki and Raymo, 2005; Clark et al., 2009); see Table 3.

When plotted in Fig. 5c against the variations in global sea level deduced from North Atlantic and equatorial Pacific δ¹⁸O benthic records by Waelbroeck et al. (2002), the length of the GIs appears to be related to the mean sea level; see Fig. S3. The long GIs occurred between 120 and 80 ka b2k and between 59 and 40 ka b2k. During these two intervals, the global sea level exhibited relatively slight variations between about −15 and −45 m and about −50 and −75 m, respectively. Conversely, after 80 and after 32 ka b2k, GIs were shorter and occurred during the most abrupt drops in sea level of the last climate cycle, from −15 to −85 m and from −75 m down to a minimum of about −120 m during the Last Glacial Maximum. The agreement between the results of the recurrence analysis and the NGRIP δ¹⁸O transitions, as well as the link between the length of GIs and the global sea level variations, seems to establish a relation between the variation in GI length and the spatial extent and elevation of the largest continental ice sheets, especially in the Northern Hemisphere, with the latter expanding further inland and out onto the continental shelf during the last climate cycle.

The DO cycles have been grouped together following an overall cooling trend of the GSs and subsequent to a strong DO event (Bond et al., 1992; Alley, 1998; Alley et al., 1999; Clark et al., 2007). This cooling trend ends with a final and coldest GS that coincides with an HE. These groupings have been named Bond cycles (Broecker, 1994; Alley, 1998), and they have mainly been observed in North Atlantic marine records of the last climate cycle with a rough periodicity of 7 kyr (Clark et al., 2007); see Fig. 6a, b. These Bond cycles, like the DO cycles, have no relationship with the periodicity of the orbital parameters (Table S1), but they are in excellent agreement with the robust 6–7 kyr periodicity of a natural, intrinsic paleoclimate oscillator (Källèn et al., 1979; Ghil and Le Treut, 1981).

Figure 6 Description of a Bond cycle. (a) Idealized Bond cycle as illustrated following Alley (1998). DO stands for Dansgaard–Oeschger event, H event stands for Heinrich event, and NADW stands for North Atlantic Deep Water. (b) Variations in the percentage of Neogloboquadrina pachyderma (s.), a species indicative of cold surface water, from DSDP 609 (Bond et al., 1992), illustrating two Bond cycles. These cycles show a series of Dansgaard–Oeschger (DO) cycles composed of an abrupt warming that is followed by a return to glacial conditions represented by “stadials”. Every Bond cycle corresponds to a long-term cooling trend that starts with a strong warming and ends with a stadial that includes a massive iceberg discharge into the North Atlantic; Heinrich events are marked by a letter “H” followed by a number assigned by Bond et al. (1992).

This oscillator is based on the countervailing effects of the positive ice–albedo feedback on Earth's radiation balance – with temperatures drops that are enhanced by the increasing extent of sea ice (Budyko, 1969; Sellers, 1969) – and of the negative precipitation–temperature feedback on the mass balance of ice sheets, with temperature increases that contribute to increased accumulation of the ice (Källèn et al., 1979; Miller and De Vernal, 1992; Ghil, 1994; Tziperman and Gildor, 2003). Ghil and Tavantzis (1983) showed that the oscillator's 6–7 kyr periodicity is quite stable over a substantial range of parameter values, and Ghil (1994) noted that this periodicity was predicted by Källén et al. (1979), well before the HEs were discovered by Heinrich (1988). More recently, the HEs' approximate recurrence time was found to be roughly equal to this periodicity by Clark et al. (2007, and references therein), in spite of the fact that these authors were not aware of its theoretical prediction by Michael Ghil​​​​​​​ and colleagues (Källèn et al., 1979; Ghil and Le Treut, 1981; Ghil, 1994).

The HEs are believed to first occur at about 0.65 Ma, as can be deduced from the U1308 bulk carbonate δ¹⁸O records. The GSs' duration has sometimes been incorrectly linked with the occurrence of an HE, leading to the misinterpretation of HEs as being equivalent to GSs. A detailed study of a subset of GSs has demonstrated that HEs did not last the entire duration of a GS (Guillevic et al., 2014), indicating much more complex dynamics in the cold stadials themselves than had initially been considered. Such complex climate behavior, by extension, may have prevailed since the first occurrence of HEs in the North Atlantic at about 0.65 Ma. Moreover, Bond and Lotti (1995) demonstrated that although an HE was embedded into the final GS of the Bond cycles, additional IRD events of lower magnitude than an HE were also embedded in the previous and intermediary GSs.

Figure 7 Schematic diagram of the proposed updated Bond cycle. (a) The amended scheme here differs from the one in Fig. 6a by adding the IRD event associated to every Greenland stadials (GS). The canonical DO are labeled here Greenland interstadials (GI). (b) The revised cycle here differs from the one in Fig. 6a by using the NGRIP δ¹⁸O data (Rasmussen et al., 2014) to mark the GSs and GIs. IRD​​​​​​​ events observed in contemporaneous marine records by Bond and Lotti (1995), are indicated by the letters “f” to “h”, while IRD events that were observed but not assigned a number also by Bond and Lotti (1995), are indicated by a letter “x”. HE numbers are the same as in Fig. 5. ​​​​​​​(c) Maps illustrating the climate evolution associated with the “long-term cooling trend” that corresponds to a Bond cycle in panel (a). The last DO cycle, named a Heinrich stadial, is characterized by a massive release of icebergs. Annual mean sea surface temperature (^∘C) for a GI (here 47 ka), a GS (here 44.4 ka), and a Heinrich stadial (at 48 ka), as simulated in a transient experiment of MIS3 (Menviel et al., 2014, 2021). GIS stands for Greenland ice sheet, LIS stands for Laurentide ice sheet, BIS stands for British Isles ice sheet, and FIS stands for Fennoscandian ice sheet.

Therefore, the previous definition of the Bond cycles given by Lehman (1993), Broecker (1994), and Alley (1998) could be revisited by including the IRD events embedded in every stadial (Fig. 7a, b). Thus, Bond cycles should be interpreted as a sequence of DO cycles that starts with a distinctly warm GI that is followed by a cooling trend with increasingly colder GSs that include as many IRD events as the stadials being identified, the latest of which is an HE; see Fig. 6a, b.

The synthetic Greenland δ¹⁸O record reconstructed from the EPICA data over the last 800 kyr by applying the bipolar seesaw model (Barker et al., 2011) indicates that DO cycles have occurred at least during this time interval. Furthermore, this synthetic record is well correlated with δ¹⁸O variations observed from Chinese speleothems (Cheng et al., 2016). Hence, the millennial variability associated with the DO cycles is likely to have existed since 0.8 Ma, or even since 0.9 Ma, when the global ice volume strongly increased, as indicated by MIS22 δ¹⁸O and sea level values. Birner et al. (2016) report millennial variability already during the much older time interval 1235–1220 ka, i.e., MIS41–37, from a marine record on the Iberian Margin. Although detecting variations in planktonic δ¹⁸O that are comparable to the MIS3 DO events in intensity and sawtooth shape, Birner et al. (2016) indicate that

identifying further Bond-like cycles in MIS38 and 40 is ambiguous. Although the lack of additional cycles might be due to the short duration of glacials in the 41 ka world, the occurrence of Bond-like cycles in the Early Pleistocene would not necessarily be expected, owing to their intrinsic relationship to Heinrich events (Bond et al., 1993) that have not been observed in the Early Pleistocene (Hodell et al., 2008).

Specifically, the closing stadial of these cycles does not show a massive IRD discharge or HE as described during the last climate cycle.

However, since IRD delivery to the North Atlantic requires ice sheets to reach the ocean, and since the first IRDs are recorded in the North Atlantic at about 1.5 Ma (Hodell and Channell, 2016a), one could assume that the start of this type of millennial variability occurred as early as this older threshold. Whether a younger start date of 0.9 Ma or an older one of 1.5 Ma is posited, these dates show that the Northern Hemisphere ice sheets played a significant role in the onset of millennial and sub-millennial climate variability that prevailed during the Middle Pleistocene and Late Pleistocene.

Ziemen et al. (2019) simulated HEs following the binge–purge model of MacAyeal (1993) and reported a two-step mechanism. The first step is a surge phase with enhanced fresh water discharge weakening the deep water formation due to the stratification of the surface water, increased sea ice cover, and leading to reduced North Atlantic sea surface temperature due to a weakened AMOC, reduced evaporation, and precipitation. The second step corresponds to a post-surge phase with a much-lowered elevation of the Laurentide ice sheet that may have lost several hundreds of meters or more during the massive iceberg discharges. During the post-surge phase described by Ziemen et al. (2019), a higher sea level associated with the lower elevation of the ice sheet favors the northern polar jet moving northwards, which leads to more precipitation over Hudson Bay, speeds up the regrowth of the Laurentide ice sheet, and results in the start of a new Bond cycle. A similar mechanism, albeit of lower magnitude, could be considered for the other iceberg discharges occurring during the GS, which did not yield HEs. Such a mechanism may apply, with reduced amplitude, to the smaller Eurasian ice sheets that were also involved in the release of icebergs during HEs. The Ziemen et al. (2019) mechanism appears, moreover, to be in agreement with the Mg $/$ Ca data on benthic foraminifera studied by Marcott et al. (2011). These authors found that warming of the subsurface temperature of the high-latitude North Atlantic led to increased basal melting under ice shelves, accelerated their collapse as suggested by Boers et al. (2018), and drove the Hudson Strait Ice Stream of the Laurentide ice sheet to release IRD deposited as HEs (Alvarez-Solas and Ramstein, 2011); see Fig. 7c.

Guillevic et al. (2014) studied ¹⁷O excess variations in Greenland ice cores, which characterize changes in the lower-latitude hydrological cycle, and reported that, at least for HE4 and HE5, iceberg delivery to the North Atlantic did not last the whole GS duration. Instead, this delivery occurred about a hundred years after the start of the GS. This short interval seems to correspond to the time during which the subsurface warming of the ocean was linked to an expansion of the sea ice and the ice shelves, as initially mentioned by Marcott et al. (2011) and Boers et al. (2018), and it also corresponds to “pre-surge” conditions. According to the Ziemen et al. (2019) model, increasing the flow of ice beyond a certain threshold leads to the massive calving of icebergs (surge phase). Furthermore, Guillevic et al. (2014) indicate that the HE4 post-surge started during GS9, hundreds of years prior to the start of the warming corresponding to GI8.

These sequences – starting with a strong GI, followed by intermediary GSs that include IRD events, and ending with a colder GS that includes an HE (see Fig. 7a,b) must have repeated throughout the last climate cycle, most likely during the glacial interval from 116 ka b2k covered in Bond and Lotti (1995). In fact, HEs are believed to have first occurred at about 0.65 Ma, although DO events may have occurred independently since about 0.8 or even 0.9 Ma. It thus appears that the millennial and sub-millennial variability corresponding to the Bond cycles involves variations in size and elevation of the Northern Hemisphere ice sheets that give rise to the iceberg release events into the North Atlantic; see Fig. 7c. This mechanism may have played a key role during the past 0.8–0.9 Myr, when the Northern Hemisphere ice sheets were at their maximum size and extended out over the continental shelf. This combination of size and contact with the much warmer ocean is likely to have destabilized the ice sheets around the North Atlantic and have led to massive IRD events. The appearance of IRD in the North Atlantic Ocean, however, might have occurred as early as 1.5 Ma (Hodell and Channell, 2016a), a fact that might suggest the prevalence of the variability over short periods discussed herein over the last 1.5 Myr.

In any case, a time interval of the most recent 0.8 or 1.5 Myr seems to have witnessed millennial-scale climate variability with an amplitude that exceeded the one due to the direct forcing by the orbital periodicities (Ghil and Childress, 1987; Ghil, 1994; Riechers et al., 2021). This fairly agreed-upon fact lends support to the interpretation of the enhanced millennial variability during glacial times as arising from an internal oscillation of the climate system – which could well be consistent with the amended Bond cycle that we have proposed – as well as with the mechanisms proposed by several authors, such as Källèn et al. (1979), Le Treut et al. (1988), Saltzman (2002) and Crucifix (2012) – while Hodell and Channell (2016a) only considered it as noise superimposed on the orbital variability deduced from their wavelet analysis of the benthic δ¹⁸O record of U1308. Substantial, nonlinear interactions between internal oscillatory variability and orbital forcing are actively being explored by Riechers et al. (2021) in this special issue.

6 Concluding remarks

Our quick overview of millennial-scale climate variability over the last 3.2 Myr suggests the following conclusions.

• The key phenomena that characterize this millennial-scale variability are Heinrich events (HEs), Dansgaard–Oeschger (DO) events, and Bond cycles. Abrupt changes are intimately interwoven with these phenomena.

• Present investigations – including both recurrence plot (RP) analysis and Kolmogorov–Smirnov (KS) methodology – point to internal mechanisms being responsible for these millennial-scale events and for the associated abrupt changes. These mechanisms include internal oscillations of the ice sheet–ocean–atmosphere system and episodic calving of ice sheets.

• The Bond cycles are linked to the dynamics of the Northern Hemisphere ice sheets, specifically to variations in their spatial extent and their elevation. The classical Bond cycles end with massive iceberg discharges into the North Atlantic Ocean mainly from the Laurentide ice sheet but also from the Fennoscandian, Greenland, Iceland, and British ice sheets. IRD releases to the North Atlantic have also been documented during every stadial. Therefore, a link with the Northern Hemisphere ice sheet extents appears evident in order to allow iceberg calving into the North Atlantic, whatever the magnitude of the calving event, i.e., either IRD events or HEs. These Bond cycles in their new interpretation illustrate a much more complex millennial variability than initially contemplated by not considering the DO and HE individually but instead linking them into a unified story.

• Millennial-scale variability is observed in proxy records from the very beginning of the last glacial period and during previous glacial periods, at least since 0.8–0.9 Ma, when HEs first appear in the records. This timing seems to coincide roughly with the MPT that has been associated with a global increase in ice volume on Earth. This coincidence does not exclude the possibility of an even earlier appearance of millennial variability (Birner et al., 2016), given the first appearance of IRD events as early as 1.5 Ma.

• Dynamical interactions between the ocean, the cryosphere with its continental ice sheets and sea ice cover, and the atmosphere are at play in generating the millennial-scale variability during glacials that leads to abrupt climate changes. However, the specific mechanisms of these interactions are still being elucidated.

• Even so, we have seen that orbital forcing, as postulated by Milankovitch, sets the stage for these internal processes and modulates their period and amplitude.

Orbital-scale and millennial variability appear to interact during the Quaternary, with millennial variability increasing in intensity, as ice sheets grew larger in both area and volume. Thus, abrupt climate changes are undoubtedly related, albeit indirectly, to the astronomical theory of climate.