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

. Abrupt climate changes are deﬁned 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, Dansgaard–Oeschger events, identiﬁed the last cycle. DO

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 14 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 δ 18 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 (1977Berger ( , 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., 2000Cheng 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., 1992Bond et al., , 1993Andrews 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.

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 . This Icehouse state is characterized by a change of the interplay between benthic δ 13 C and δ 18 O, which corresponds to a new relationship between the carbon cycle and climate. Indeed, the negative excursions in δ 13 C were associated with negative ones in δ 18 O during most of the Cenozoic since 66 Ma. A shift occurred, however, in the Plio-Pleistocene at about 5 Ma, with negative δ 13 C excursions being associated with positive δ 18 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.
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 δ 18 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 δ 18 O record best describes the more or less regular recurrence of cold and warm events that we are analyzing; see Fig. 1c.

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 δ 18 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. (2007Marwan et al. ( , 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 paleorecords. 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 δ 18 O or the GS-GI boundaries in the North Greenland Ice Core Project (Barbante et al., 2006) 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. (2007Marwan et al. ( , 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 = 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 https://doi.org/10.5194/cp-18-249-2022 Clim. Past, 18, 249-271, 2022 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).
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 δ 18 O, 0.127 for the U1308 bulk carbonate δ 18 O, and 0.173 for NGRIP δ 18 O. The results were plotted along with the original recurrence plot (Table 1). Please see further details in Bagniewski et al. (2021).

The past 3.2 Myr history of the Northern
Hemisphere climate Analyzing the U1308 benthic δ 18 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 δ 18 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 Qua-ternary, 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 δ 18 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 10 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 δ 18 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), 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).
The benthic δ 18 O record of the U1308 marine-sediment core is interpreted in terms of global ice volume and deepocean temperatures (Chappell and Shackleton, 1986;Elderfield et al., 2012). Its recurrence analysis shows a drift topology (Marwan et al., 2007(Marwan et al., , 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  (Hodell and Channell, 2016a), and NGRIP δ 18 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.  (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. The δ 18 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 δ 18 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 δ 18 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 . The reconstructed CO 2 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 . After 1.25 Ma, the sea level changes increased to about 70-120 m below the present-day value, while the reconstructed CO 2 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 2 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.

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).
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   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;Menviel et al., 2014;Ditlevsen, 2018, 2019) it corre-  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/. sponds 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., 2000Müller et al., 2003;Fletcher et al., 2010;Rousseau et al., 2017a, Table 3. Thresholds identified in the recurrence plot of NGRIP δ 18 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 (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 δ 18 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) 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 δ 18 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.

260
D.-D. Rousseau et al.: Abrupt climate changes and the astronomical theory: are they related?

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 15 N measurements from the Greenland ice . 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., 2017aRousseau et al., , 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 , 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. The recurrence analysis performed on the high-resolution NGRIP δ 18 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 δ 18 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 δ 18 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 δ 18 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 δ 18 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).
This oscillator is based on the countervailing effects of the positive ice-albedo feedback on Earth's radiation balancewith 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 δ 18 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 , 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.
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 (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).
The synthetic Greenland δ 18 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 δ 18 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 δ 18 O and sea level values. Birner et al. (2016) report millennial variability already dur-ing 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 δ 18 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  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 bingepurge 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  Fig. 6a by using the NGRIP δ 18 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(Menviel et al., , 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. 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 17 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 agreedupon 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 δ 18 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.

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 millennialscale 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.
Author contributions. DDR designed the project and drafted the paper. WB performed the recurrence analyses. All of the authors contributed to the writing of the paper. DDR and MG finalized the paper.
Competing interests. At least one of the (co-)authors is a member of the editorial board of Climate of the Past. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.
Disclaimer. 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 "A century of Milankovic's theory of climate changes: achievements and challenges (NPG/CP inter-journal SI)". It is a result of the conference "One Hundred Years of Milankovic's Theory of Climate Changes: synergy of the achievements and challenges of the next century", 17-18 November 2020.
Acknowledgements. The authors thank the editor, two anonymous reviewers, and Linda Hinnov for their handling of the paper and useful comments that greatly benefited the revised version of the original draft. It is a pleasure to acknowledge Laurie Menviel for providing the sea surface maps used in Fig. 6 and Peter Clark for his comments on the draft of this paper. An earlier version of this paper was presented at the Milankovitch symposium held in 2020 to celebrate the centennial of the Milankovitch (1920)