Abrupt cl imate changes and the astronomical theory

Abrupt cl imate changes and the astronomical theory Denis-Didier Rousseau, Witold Bagniewski, and Michael Ghil 1 Geosciences Montpellier, University Montpellier, CNRS, Montpellier, France 2 Lamont-Doherty Earth Observatory, Columbia University, New York, 10964, USA 3 Laboratoire de Météorologie Dynamique (Institut Pierre Simon Laplace and CNRS), Ecole Normale 5 Supérieure and PSL University, Paris, France 4 Department of Atmospheric and Oceanic Sciences, University of California, Los Angeles, CA 90095, USA


Introduction
Well-dated geological data indicate that the Earth experienced orbitally paced climate changes since at least the late Precambrian -1.4 billion years ago during the Proterozoic Eon (Benn et al., 2015;Zhang et al., 2015; 40 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., 2017Kent et al., , 2018;;Olsen et al., 2019;Drury et al., 2020;Westerhold et al., 2020).These changes reflect the variations in the Earth's axis of rotation -precession and tilt, and in the geometry of the Earth's orbit around the sun -eccentricity, driven by gravitational interactions within the Solar system (Berger, 1977;Laskar et al., 2011).These variations affect the distribution of insolation 45 at the top of the atmosphere, forcing latitudinal and seasonal climate changes with periodicities of tens or hundreds of thousands of years.Although well acknowledged presently in the climate community, this mechanism took a long journey to reach such general acceptance (Imbrie & Imbrie, 1986).
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).

50
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 keep temperatures cold enough 55 to allow ice sheets to build up.
Five decades later, Croll (1890) presented his theory of the ice ages being driven by the changing distance between the Earth and the Sun as measured on December 21, 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 60 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 65 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.

70
Thirty years after Croll, Milankovitch (1920, 1941), inspired by his exchanges with W. Köppen and A. 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 today his name, implies a fairly regular, multi-thousand-year variability.Based on this crucial idea and on using the recent orbital calculations 75 of Pilgrim (1904), among others, Milankovitch was able to estimate temperatures and insolation at various latitudes at the top of the atmosphere and to 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 Penck and Brückner's (1909) study of Alpine glaciations.
Milankovitch's astronomical theory of climate was severely criticized by contemporary physicists and rejected 80 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),

85
linked more reliably 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 Islands 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 Ma up to 1.25 Ma, and the dominant 100-kyr cycles of the Mid-and Late, or Upper, Pleistocene, i.e. the last 800 kyr that show longer glacials implying much larger 90 continental ice sheets.The transition interval, between 1.25 Ma 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 still a matter of debate (Ghil & Childress, 1987, Sec. 12.7;Ghil, 1994;Saltzman, 2002;Clark et al., 2021).

95
Although the broad astronomic framework for past climate changes seems to be widely accepted, recent highresolution investigations in ice, marine and terrestrial records revealed much shorter periodicities than the orbital ones, as well as abrupt changes that do not match those attributed to variations in eccentricity, obliquity, or precession of the equinoxes.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 100 framework is shown herein to be still quite relevant to the study of abrupt and large climate changes.
The paper is organized as follows.In Sec. 2, proxy records for the last 3.2 Myr of Northern Hemisphere climate are described.In Sec. 3, we concentrate on the millennial-scale variability revealed by these records.In Sec. 4, we establish a connection between Dansgaard-Oeschger (DO) events and amended  (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 δ 13 C and δ climate are particularly well described in North Atlantic core U1308 (Hodell and Channell, 2016).This core is located in the ice-rafted debris (IRD) belt (Ruddiman, 1977), and it yields a more complete record than U1313 (Naafs et al., 2013); see Fig. 125 variations, and they occur at 2.75 Ma, 1.5 Ma, 0.9 Ma 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 in the ocean, and the melting of these icebergs is likely to have impacted the oceanic circulation.Naafs et al. (2013), however, reported the occurrence of weak IRD events in the late Pliocene that they attributed mainly to 130 Greenland and Fennoscandian glaciers, and that point to larger ice sheets over these regions than during the later Quaternary, when North American ice sheets were considerably larger.
The second date corresponds to an increased amplitude in ice volume variations between glacial minima and interglacial optima.This second step shows the permanent occurrence of ice-rafted events during glacial intervals in the record, therefore an amplified relationship of climate variations with Northern Hemisphere ice 135 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 Antractic 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 140 the Alps (Knudsen et al., 2020;Muttoni et al., 2003).At about the same time, the synthetic Greenland δ 18 O reconstruction indicates the occurrence of millennial variability expressed by DO-like events (Barker et al., 2011).
The last date at 0.65 Ma marks the end of the transition from the Lower and Mid-Pleistocene intervalcharacterized by 41-kyr-dominated cycles and smaller 23-kyr ones -to the Upper Pleistocene, with its 100-145 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, especially due to the larger contribution of the Northern American ice sheets.The latter now have a bigger impact on Northern 150 Hemisphere climate than the Eurasian ice sheets (Batchelor et al., 2019).The IRD event intensity and frequency of occurrence increase (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 Ma -0.4 Ma is also the interval during which Northern Hemisphere ice sheets reached their southernmost extent (Batchelor et al., 2019).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 = 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 !} of the series on the abscissa and another copy {x !} on the ordinate.A dot is entered into a position (i, j) of the matrix when x ! is sufficiently 165 close to x ! .For the 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, in general, the matrix is rather symmetric, although one does not always define closeness symmetrically; to wit, x ! may be "closer to" x !than x ! is to x !(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 170 latter is certainly the case for the climate system on time scales of 10-100 kyr and longer, which is 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 175 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 exactly periodic but only nearly so.Hence, such patterns are not that 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 more objectively quantifiable by recurrence quantification analysis and bootstrapping (Efron, 1981; 180 Efron and Tibshirani, 1986).
The benthic δ 18 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) associated with nonstationary systems with slowly varying parameters.The RP exhibits, moreover, a characteristic texture -given by the pattern of  (Shackleton and Opdyke, 1977;Pisias and Moore, 1981;Ruddiman et al., 1989;Clark and Pollard, 1998;Clark et al., 2006).Channel, 2016).The recurrence analysis of this record also displays a drift topology, and it yields the two 195 Hodell & Channel (2016) steps at 1.5 Ma and 0.65 Ma (Fig. 3a, Tab. 1).Our analysis further identifies the steps at 0.9 Ma and 2.75 Ma, however, it also detects other steps at 1.25 Ma (also noticed for the δ 18 O).Table 1, comparing Hodell & Channel (2016) 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-

215
The behavior of the U1308 proxy records on the time scale 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.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. 220 3).
Observations of such abrupt variations have been reported in some detail 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., 225 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 230 does not exceed a millennial time scale (Broecker, 1994;Boers et al., 2018;Boers, 2018).As the periodicity of the events is at millennial and submillenial scale (Broecker, 1994;Clark et al., 1999;Ganopolski and Rahmstorf, 2001;Rahmstorf, 2002;Schulz, 2002;Menviel et al., 2014;Lohmann andDitlevsen, 2018, 2019) 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.

235
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).Dansgaard et al (1993), though, clearly identified 23 rapid 240 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), and they are considered as the "canonical" DOs.
The 23 canonical DOs were described later using various types of record, in marine sediments (Bond et al., 245 1992;Henry et al., 2016), as well as in terrestrial ones (Allen et al., 1999;Sanchez-Goni et al., 2000, 2002;Müller et al., 2003;Fletcher et al., 2010;Rousseau et al., 2017Rousseau et al., a, b, 2021)), including speleothems (Wang et al., 2001;Genty et al., 2003;Fleitmann et al., 2009;Boch et al., 2011).The increase in resolution of the Greenland ice core investigations (Fischer et al., 2015;Schupbach et al., 2018;Svensson et al., 2020) and of several speleothems (Cheng et al., 2016) allowed one to conduct much more detailed analyses of the δ 18 O but also of 250 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 GIs, followed by their associated GSs.
In addition to the unresolved dispute on defining this interstadial and stadial labeling (Rousseau et al., 2006), remains the status of and significance attributed to the smaller-scale events detected by Rasmussen et al. (2014).Bagniewski et al. (2020Bagniewski et al. ( , 2021) ) recently proposed a method based on an augmented nonparametric and 87 ka b2k, had on average a range of 10°C -12°C (Kindler et al., 2014), with each transition lasting on average between 50 yr and 100 yr (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 (Boers et al., 2018),

295
but the timing does not correspond to any periodicity of the orbital parameters, neither in the average duration of the events, nor in the average interval between two such events.
The recurrence analysis performed on the high-resolution NGRIP δ 18 O record and illustrated in Fig. 4b 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 Fig. 3, the NGRIP study shows a drift topology but with a non-300 uniform pattern.Several thresholds are identified which 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 Tab.The DO cycles have been grouped together following an increasing cooling trend of the GSs and subsequent to 315 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 been observed mainly in North Atlantic marine records of the last climate cycle with a rough periodicity of 7 kyr (Clark et al., 2007); see Fig. 5.These Bond cycles, like the DO cycles, have no relationship with the periodicity of the orbital parameters, but they are in excellent agreement 320 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 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,

325
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 equal roughly to this 330 periodicity by Clark et al. (2007, and references therein), in spite of the fact that these authors were not aware of its theoretical prediction of M. Ghil and coauthors (Källén et al., 1979;Ghil and Le Treut, 1981;Ghil, 1994).
The HEs events are believed to first occur at about 0.65 Ma, as can be deduced from the U1308 benthic and 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 335 has demonstrated that HEs did not last the entire duration of a GS (Guillevic et al., 2014), indicating a much more complex dynamics of 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 340 previous and intermediary GSs.
Therefore, the previous definition of the Bond cycles given by Broecker (1994) and Alley (1998), could be amended by including the IRD events embedded in every stadial.Thus Bond cycles should be interpreted as a sequence of DO cycles, which start with a distinctly warm GI, followed by a cooling trend with increasingly colder GSs that include as many IRD events as identified stadials, the latest of which is an HE; see Fig. 6a.

345
The synthetic Greenland δ 18 O record reconstructed from the EPICA data over the last 800 ka by applying the bipolar seesaw model (Barker et al., 2011) indicates that DO cycles have occurred at least during this time interval.This synthetic record is, furthermore, 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 prevailed since 0.8 Ma, or even since 0.9 Ma, when the global ice volume strongly increased, as indicated Channell, 2016), 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.5Ma, is posited, these results show that the Northern Hemisphere ice sheets played a significant role in the onset of millennial and sub-

355
millennial climate variability that prevailed during the Mid and Late Pleistocene.

375
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 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 380 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. 6a -must have repeated throughout the last climate cycle, most likely so during the glacial interval from 116 ka b2k on (Bond and Lotti, 1995).In fact, HEs are 385 believed to first occur at about 0.65 Ma, although DO events may have occurred independently since about 0.8 Ma 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. 6b.This mechanism may have played a key role during the past 0.8 Myr-0.9Myr, when the Northern Hemisphere ice sheets were at their maximum size and extended

395
In any case, a time interval of the most recent 0.8 Myr-1.5 Myr seems to have witnessed millennial-scale climate variability whose amplitude exceeded the one that was directly forced by the orbital periodicities (Ghil and Childress, 1987;Ghil, 1994;Ghil, 2021).This fairly well-agreed-upon fact leads support to the interpretation of the enhanced millennial variability during glacial times as arising from an internal oscillation of the climate system -as proposed by several authors, such as Källèn et al. (1979), Le Treut et al. (1988), Saltzman (2002) 400 and Crucifix (2012), among others -while Hodell and Channell (2016) 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) 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 millennial-scale variability are Heinrich events (HEs), Dansgaard-Oeschger (DO) events, and Bond cycles.Abrupt changes are intimately interwoven with 410 these phenomena.
• Present investigations 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 415 variations in their spatial extent and their elevation.These Bond cycles illustrate a much more complex millennial variability than initially contemplated, by considering the DO and HE not individually, but 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.9Ma, when HEs first appear in the 420 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, 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 that leads to abrupt 425 climate changes.The specific mechanisms of these interactions, though, 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

Datasets used
The datasets used include the following:       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) Stadial, is characterized by a massive release of icebergs.Annual mean SST (°C) for a GI (here 47 ka), a GS (here 44.4 ka), and a Heinrich stadial (here 48 ka), as simulated in a transient experiment of MIS3 (Menviel et al., 2014(Menviel et al., , 2021)).GIS = Greenland Ice Sheet, LIS = Laurentide Ice Sheet, BIS = British Isles Ice Sheet, and FIS = Fenno-Scandian Ice Sheet. 940 https://doi.org/10.5194/cp-2021-103Preprint.Discussion started: 11 August 2021 c Author(s) 2021.CC BY 4.0 License.
155 https://doi.org/10.5194/cp-2021-103Preprint.Discussion started: 11 August 2021 c Author(s) 2021.CC BY 4.0 License.The descriptions and classifications of climatic phenomena discussed so far were essentially subjective, being based on the mere visual inspection of the proxy records and on the previous experience of the investigators with the study of similar records.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 the recurrence plots (RPs) introduced by Eckmann et al. (1987) into the study of dynamical systems and popularized in the climate 160

185
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 δ 18 O variability (Fig.2a, Tab. 1).Two are roughly similar to those determined by Hoddell and Channel (2016): at 1.5 Ma and 0.65 Ma, and three differ: at 2.95 Ma, 2.55 Ma, and 1.25 Ma.Interestingly, the interval 1.25 Ma 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 190 100-kyr dominated ones occurred The δ 18 O bulk carbonate record in the U1308 core, in turn, is interpreted as characterizing IRD released into the North Atlantic Ocean, with the most negative δ 18 O values representing the largest iceberg calvings (Hodell & https://doi.org/10.5194/cp-2021-103Preprint.Discussion started: 11 August 2021 c Author(s) 2021.CC BY 4.0 License.

200
interglacial fluctuations.The interval 2.8 to 1.2 Ma shows glacial-interglacial sea level variations of about 25-50 m below the present day.The CO 2 concentrations varied between 270 ppmv and 280 ppmv during interglacials and between 210 ppmv 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 decreased to about 70-120 m below the present day, while the CO 2 concentrations varied between 250 ppmv and 320 ppmv during 205 interglacials and between 170 ppmv 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 "Milanković glacials," which correspond to the odd marine isotope stages determined in the U1308 core and in many others, are characterized by low eccentricity and obliquity, and a boreal summer that coincides with aphelion and leads, 210 therefore, to minimum values of summer insolation.The increase in IRD variability and magnitude since 1.5 Ma, however, shows that distinct, faster processes have to be considered than those due to slow changes in Earth's orbital parameters; see again Figs. 2 and 3. 3 Millennial-Scale Variability , it https://doi.org/10.5194/cp-2021-103Preprint.Discussion started: 11 August 2021 c Author(s) 2021.CC BY 4.0 License.

255
Kolmogorov-Smirnov (KS) test to detect abrupt transitions from paleo-records in a highly robust way.The classical KS test consists of quantifying and comparing the empirical distribution functions of two samples from a time series, before and after a potential jump.In the present case, the test is augmented by varying window sizes, and evaluating the rate of change and the trend in maxima and minima of the time series, thus establishing the main transitions in a record, such as the GS-GI boundaries in the North Greenland Ice Core 260 Project (NGRIP) δ 18 O.The transitions identified by this methodology include all the canonical events described by Rasmussen et al. (2014), although a few other DOs or subevents are not detected by it, even after changing the test window size.Using a global climate indicator like methane (CH 4 : Blunier and Brook, 2001), EPICA Community Members (2006) and the WAIS Consortium (2015) demonstrated that the millennial-scale variations observed during the 265 last 130 kyr in Greenland are observed in Antarctica as well, and that they are, therefore, a global phenomenon.The δ 18 O variations in two hemispheres, though, are in opposite phases, with Southern Hemisphere warmings occurring prior 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 270 Circulation (AMOC) play a key role in this teleconnection as AMOC slowdown during GSs corresponds to reduced northward heat transport (Sarnthein et al., 2001; Ganopolski and Rahmstorf, 2001; McManus et al., 2004).Detailed high-resolution marine-sediment and ice core studies (Buizert et al., 2015a; Henry et al., https://doi.org/10.5194/cp-2021-103Preprint.Discussion started: 11 August 2021 c Author(s) 2021.CC BY 4.0 License.2016), on the other hand, indicate a southward propagation of the anomalies, with Greenland climate leading Antarctica by approximately 200 years.At this point, one can only state that the latter, 275 observational studies do not contradict the earlier, modeling studies, but do not directly confirm them either.To address this issue, Boers et al. (2018) recently developed a simple model to reconstruct the millennial variability in δ 18 O of the past 60 kyr b2K, as observed in the high-resolution ice cores from NGRIP in Greenland and West Antarctic Ice Sheet (WAIS) in Antarctica.This simple model, based on the bipolar seesaw 280 mechanism of 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, the 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, whether subsumed by the 65°N summer 285 insolation curve or not.4. 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 (Guillevic 290 et al., 2014).These estimates indicate that the 21 warming events during the last climate cycle, between 12 ka 2. When plotted in Fig. 4c 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 305 mean sea level.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 m and -45 m and about -50 m and -75 m, respectively.Conversely, after 80 ka b2K 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 m to -85 m, and from -50 m down to a minimum of about -120 m during the Last Glacial Maximum.The agreement between the 310 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 https://doi.org/10.5194/cp-2021-103Preprint.Discussion started: 11 August 2021 c Author(s) 2021.CC BY 4.0 License.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.
350 by MIS 22 δ 18 O and sea level values.However, since IRD delivery to the North Atlantic requires ice sheets to reach the ocean, and since the first IRD are recorded in the North Atlantic at about 1.5 Ma (Hodell and https://doi.org/10.5194/cp-2021-103Preprint.Discussion started: 11 August 2021 c Author(s) 2021.CC BY 4.0 License.
Ziemen et al. (2019) simulated HEs following the Binge-Purge model ofMacAyeal (1993) and reported a twostep 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 360second 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 postsurge phase described byZiemen et al. (2019), a higher sea level associated with the lower elevation of the ice sheet favors the northern polar jet to move 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 365 mechanism, albeit of lower magnitude, could be considered for the other iceberg discharges occurring during GS, which did not yield HEs.Such a mechanism may apply, with a 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 byMarcott et al. (2011).These authors found that warming of the subsurface temperature of the high-latitude North Atlantic lead 370 to increased basal melting under ice shelves, accelerated their collapse as suggested byBoers et al. (2018), and thus drove the Hudson Strait Ice Stream of the Laurentide ice sheet to release IRD deposited as HEs (Alvarez-Solas and Ramstein, 2011); see Fig.6b.Guillevic et al. (2014) studied17 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 390 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 led to massive IRD events.The appearance of https://doi.org/10.5194/cp-2021-103Preprint.Discussion started: 11 August 2021 c Author(s) 2021.CC BY 4.0 License.IRD in the North Atlantic Ocean, though, might have occurred as early as 1.5 Ma (Hodell & Channel, 2016), a fact that could indicate the type of relatively short-periodic variability discussed herein prevailing over the last 1.5 Myr.
430 undoubtedly related, albeit indirectly, to the astronomical theory of climate.https://doi.org/10.5194/cp-2021-103Preprint.Discussion started: 11 August 2021 c Author(s) 2021.CC BY 4.0 License.for inviting the lead author.DDR, WB and MG are funded by the European Union's Horizon 2020 research and innovation programme through the TiPES grant no.820970.This is LDEO contribution (no.ZZZ), and TiPES contribution (no.x).Financial support 450 This research has been supported by the European Union's Horizon 2020 research and innovation programme (TiPES grant no.820970).

Figure captions Figure 1 .Figure 2 .
Figure captions Figure 1.Earth Climate history represented by two marine records and an ice core one.(a) Record of benthic 895

Figure 3 .
Figure 3. Recurrence analysis of U1308 bulk carbonate δ 18 O.(a) RP of the bulk carbonate δ 18 O, with the same proximity threshold as in fig.2A.(b) same as fig.2B, except that the vertical bars now represent the seven thresholds determined by the analysis of the time series in fig.3a.The RP web site is http://www.recurrence-910

Figure 6 .
Figure 6.Schematic diagram of the proposed, updated Bond cycle.(a) The amended scheme here differs from the one in fig. 5 by using the NGRIP δ 18 O data of Rasmussen et al. (2014) to mark the Greenland stadials 930 c g x Z o A w w e w T N 4 B W / G k / F i v B s f s 9 E l o 9 g 5 A H 9 g f P 4 A o y a X T g = = < / l a t e x i t > ( 18 O) https://doi.org/10.5194/cp-2021-103Preprint.Discussion started: 11 August 2021 c Author(s) 2021.CC BY 4.0 License.
, are indicated by a letter 'x'.HE numbers as in fig. 5. (b) Maps illustrating the climate evolution associated with the "Long-Term Cooling Trend" that corresponds to a Bond cycle in panel (A) above; the last DO cycle, named a Heinrich 935

Fig. 1
Fig.1 Fig.2 Fig.5960 Bond cycles, and outline how global ice sheet volume affects the latter.Concluding remarks follow in Sec. 5.
last two succeeded each other, thus generating the classical climate trend towards the recent ice age conditions 1b.The variations of the benthic δ 18 O mostly indicate varying periodicities through time that correspond to periodicities in the orbital parameters of Earth's climate, as also pointed out by Lisiecki

Table captions Table 1 .
Comparison of the main steps detected by Hodell and Channell (2016) from U1308 marine record with those deduced from the Recurrence Plot of the benthic δ 18 O and bulk carbonate δ 18 O data of the same record.