Mode transitions in Northern Hemisphere glaciation : co-evolution of millennial and orbital variability in Quaternary climate

We present a 3.2 Myr record of stable isotopes and physical properties at IODP Site U1308 (reoccupation of DSDP Site 609) located within the ice-rafted detritus (IRD) belt of the North Atlantic. We compare the isotope and lithological proxies at Site U1308 with other North Atlantic records (e.g., sites 982, 607/U1313, and U1304) to reconstruct the history of orbital and millennial-scale climate variability during the Quaternary. The Site U1308 record documents a progressive increase in the intensity of Northern Hemisphere glacial–interglacial cycles during the late Pliocene and Quaternary, with mode transitions at ∼ 2.7, 1.5, 0.9, and 0.65 Ma. These transitions mark times of change in the growth and stability of Northern Hemisphere ice sheets. They also coincide with increases in vertical carbon isotope gradients between the intermediate and deep ocean, suggesting changes in deep carbon storage and atmospheric CO2. Orbital and millennial climate variability co-evolved during the Quaternary such that the trend towards larger and thicker ice sheets was accompanied by changes in the style, frequency, and intensity of millennial-scale variability. This coevolution may be important for explaining the observed patterns of Quaternary climate change.


Introduction
Earth's climate during the last 2.7 Myr has been has been characterized by the waxing and waning of large Northern Hemisphere continental ice sheets.Long sediment records from the North Atlantic Basin recovered by Integrated Ocean Drilling Program (IODP) and its preceding programs (Deep Sea Drilling Project, DSDP, and Ocean Drilling Program, ODP) have provided a detailed history of Northern Hemisphere glaciation during the Pliocene and Quaternary.The intensification of Northern Hemisphere glaciation began in the latest Pliocene (∼ 2.7 Ma) and the intensity, shape, and duration of glacial-interglacial cycles changed during the Quaternary.The average climate state evolved towards generally colder conditions with larger and thicker ice sheets, and the spectral character of climate variability shifted from a dominant period of 41 kyr to a quasi-period between 80 and 120 kyr.Glaciations were generally less intense, shorter in duration, and more symmetrical during the "41 kyr world" of the early Pleistocene (Maslin and Brierly, 2015).Across the Middle Pleistocene transition (MPT), the glacial cycle lengthened, ice volume increased, and the shapes of marine isotopic stages assumed a more asymmetric, saw-tooth pattern during the "100 kyr world" of the late Pleistocene.The transition is often viewed as a shift from a linear response of Earth's climate system to obliquity forcing prior to the MPT to a more nonlinear response afterwards, although the system response was likely more complicated (e.g., Ashkenazy and Tziperman, 2004).
The causes of these long-term patterns of Quaternary climate have been attributed to internal changes in climate response because orbital forcing did not change significantly over this time.Many explanations have been invoked, including the following: 2. changes in global deep-ocean circulation resulting in modification of deep-ocean carbon storage capacity and heat distribution (Hodell and Venz-Curtis, 2006); 3. glacial erosion and associated changes in ice-sheet dynamics (Pisias and Moore, 1981;Berger and Jansen, 1994;Clark and Pollard, 1998); 4. sea ice switch mechanism as a result of the gradual cooling of the deep ocean during the Pleistocene (Tziperman and Gildor, 2003); 5. ice-sheet behavior as a multi-stable dynamical system with bifurcation points (Abe-Ouchi et al., 2013, Ditlevsen, 2009); 6. stochastic variability (Meyers and Hinnov, 2010;Ditlevsen, 2009;Huybers, 2009) None of these explanations are mutually exclusive and all could contribute by varying degrees to the observed patterns of Quaternary climate change.
Although it is accepted that orbitally induced changes in insolation act as the pacemaker of the ice ages (Hays et al., 1976), we still lack a complete understanding of what caused glacial-interglacial cycles (Raymo and Huybers, 2008;Paillard, 2015).This uncertainty is possibly because the nonlinear climate system responds not only to longer-term external and internal forcing but also to events (triggers) that can result in major reorganization of the ocean-atmosphere system (Berger, 2013;Broecker and Denton, 1989).Thus, it is important to understand if and how short-term (millennial) and long-term (orbital) climate variability interact to produce the observed patterns of Quaternary climate change.
Observations of abrupt climate change in Greenland, beginning in the early 1990s (e.g., Dansgaard et al., 1993), sparked a proliferation of studies of millennial-scale climate variability for the last glacial cycle that is now being extended to older parts of the Quaternary.The leading mechanism to explain millennial-scale oscillations in the North Atlantic is freshwater forcing of the strength of thermohaline circulation (Broecker and Denton, 1989, among others), although other processes may also be involved (Barker et al., 2015).The North Atlantic is one of the most climatically sensitive regions in the world ocean because of its proximity to the North American, Greenland, and European ice sheets.Most of the water stored in Northern Hemisphere ice sheets during Quaternary glaciations was discharged into the Atlantic Ocean, either directly or indirectly via the Arctic (Fig. 1).The buildup of ice on Northern Hemisphere continents can be thought of as a "capacitor" that stores freshwater on land during times of ice growth and releases it to the ocean during times of ice decay as icebergs and meltwater (Bender, 2013).The volume, rate, and location of freshwater discharge to the North Atlantic Ocean relative to the source areas of deepwater formation can have a strong impact on Earth's global climate.As ice sheets grew larger during glacial periods of the Quaternary, the freshwater capacitor became more highly charged and the potential for strong climate response upon discharge was enhanced.The changing topography of the ice sheets also affected the wind field (Wunsch, 2006) and the coupled atmosphere-ocean system (Zhang et al., 2014).Understanding the history of changes in the volume, height, and location of ice buildup on Northern Hemisphere continents and their impact on atmospheric and ocean circulation is important for understanding Quaternary climate evolution.
Here we present a 3.2 Myr record of stable isotopes and physical properties at IODP Site U1308 (49 • 52.67 N, 24 • 14.289 W; Fig. 1), which is located within the icerafted detritus (IRD) belt of the North Atlantic (Ruddiman, 1977).Site U1308 represents the reoccupation of ODP Site 609, which has played an important role in understanding Quaternary orbital and millennial climate change, including the recognition of Heinrich events and correlation of millennial-scale climate variability between marine sediment and Greenland ice cores (Broecker et al., 1992;Bond et al., 1992Bond et al., , 1993Bond et al., , 1999;;Bond and Lotti, 1995).We integrate the isotope and lithological proxies from Site U1308 with other North Atlantic records (e.g., sites 982, 607/U1313, and U1304) to elucidate the patterns of orbital and millennialscale variability in the subpolar North Atlantic during the Quaternary.

Composite section
Six holes were drilled at Site U1308 to ensure complete recovery of the stratigraphic section, and a shipboard composite splice was constructed to 248 mcd (Expedition 303 Scientists, 2006a).We adopted a revised meters composite depth (rmcd) scale that differs slightly from that used by Hodell et al. (2008) and Channell et al. (2016;Supplement Table S1).

Chronology
The age-depth control points for Site U1308 are given in Table S1 and plotted with sedimentation rates in Fig. S1 in the Supplement.The integrated oxygen isotope and magnetostratigraphy of Site U1308 is described elsewhere for the interval from 0 to 1.5 Ma (Hodell et al., 2008;Channell et al., 2008) and from 1.5 to 3.2 Ma (Channell et al., 2016).The interval younger than 76 ka was dated using the age model of Obrochta et al. (2012Obrochta et al. ( , 2014)), which is based on correlating variations in sediment lightness between sites U1308 and 609 and then transferring the Site 609 age model to U1308.The age model of Site 609 consists of recalibrated radiocarbon dates and correlation of Neogloboquadrina pachyderma (sin) to Greenland ice-core δ 18 O placed on the GICC05 age model.For the interval from 76.5 to 110.5 ka, ages were derived by correlation of sediment lightness at Site U1308 to NGRIP on the GICC05 timescale (S.Obrochta, personal communication, 2016).Beyond 110.5 ka, the age model was derived by correlating the benthic δ 18 O signal to the LR04 benthic oxygen isotope stack (Lisiecki and Raymo, 2005) assuming linear sedimentation rates between tie points.Modification to the oxygen isotope age model of Hodell et al. (2008) was made near the transition of Marine Isotope Stage (MIS) 10 to MIS 9 (Termination IV), where benthic foraminifera are very scarce and the original age model is inaccurate.In this interval, analysis of δ 18 O from planktonic foraminifera permitted a refined stratigraphy.The oxygen isotope record of Site U1308 indicate that the section is complete except for a short hiatus that removed MIS G1 and G2 (Fig. 2), corresponding to the interval from ∼ 2.6 to 2.65 Ma (Channell et al., 2016).

Stable isotopes
Foraminifera were picked from the > 212 µm size fraction, and one to five individuals were used for analysis.Stable oxygen and carbon isotopes were measured on the benthic foraminifer Cibicidoides wuellerstorfi and/or Cibici-doides kullenbergi using the method described in Channell et al. (2016).We measured 312 pairs of C. wuellerstorfi and C. kullenbergi from the same samples at Site U1308 and found no consistent offset for either δ 18 O or δ 13 C (Fig. S2).The sample spacing for foraminifer stable isotopes was approximately every 2 cm for the upper 100 mcd (Hodell et al., 2008) and 5-10 cm for the 100 to 248 rmcd interval of Site U1308.The analytical error was estimated to be ±0.08 ‰ for δ 18 O and ±0.03 ‰ for δ 13 C by measuring eight standards (NBS-19) with each set of 38 samples.
Oxygen isotopes of bulk carbonate were measured using a Thermo Scientific GasBench II, equipped with a CTC autosampler coupled to a Delta V mass spectrometer (Spötl and Vennemann, 2003).Analytical precision is estimated to be ±0.1 ‰ for δ 18 O by repeated analysis of the Carrara marble standard.The sample spacing for bulk oxygen isotope measurements was approximately every 2 cm for the upper 100 mcd (Hodell et al., 2008) and 4 cm for the 100-248 mcd interval of Site U1308, which corresponds to an average temporal spacing of 250 and 500 years, respectively.For Site 982, bulk samples were ground and measured at a sample spacing of every 5 cm using a Finnigan MAT Kiel III carbonate preparation device attached to a Finnigan MAT 252 mass spectrometer at the University of Florida.All isotope results are reported relative to Vienna Pee Dee Belemnite (VPDB).

Physical properties
Measurements of density and natural gamma radiation (NGR) were made onboard the R/V JOIDES Resolution during IODP Expedition 303 (Expedition 303 Scientists, Methods, 2006b).NGR was measured at a sample spacing of either 2.5 or 5 cm and density at 2.5 cm spacing.Because of the response curve of the NGR detectors is ∼ 17 cm (full width at half maximum), the record is only suitable for studying orbital-scale variation.Volume susceptibility was measured at 1 cm intervals on U-channel samples (2 × 2 × 150 cm 3 continuous samples) using a susceptibility bridge designed for U-channel samples with a response function half-peak width of 3 cm (Thomas et al., 2003).The magnetic grain size parameter (κ ARM /κ) was measured every 1 cm, but each measurement is not independent of adjacent measurements owing to the ∼ 4.5 cm width at half height of the magnetometer response function (Channell et al., 2016).
The temporal resolution of measurements of physical properties varies from 250 to 625 years assuming an average sedimentation rate of 8 cm kyr −1 for the entire record.Although interval sedimentation rates vary considerably, the resolution of the physical properties should be sufficient to detect millennial events, as demonstrated for Heinrich and other IRD events in the late Pleistocene (Hodell et al., 2008;Channell et al., 2012).Lisiecki and Raymo, 2005) and selected glacial marine isotope stages are labeled.The oxygen isotope stratigraphy is complete except for a hiatus that has removed MIS G1 and G2 (∼ 2.6 to 2. 65 Ma).Inclination of natural remanent magnetization (NRM) with polarity chrons, subchrons, and excursions labeled (Channell et al., 2016).

Time series analysis
Traditional time series analysis was conducted using RED-FIT, and spectral peaks were evaluated against a red-noise background from an AR1 process (Schulz and Mudelsee, 2002).To track the time-varying amplitude of orbital and suborbital periods, we calculated the continuous wavelet transform using the MATLAB code of Grinsted et al. (2004).Time series for wavelet analysis were prepared by Gaussian interpolation to a fixed time increment of 1 kyr using a 3 kyr window.The statistical significance of wavelet power was tested relative to a red-noise background power spectrum.

Oxygen isotopes
We compare the benthic δ 18 O record of Site U1308 to Site 607 for the interval from 0 to 2. 4 Ma and to Site U1313 (re-occupation of 607) for 2.4 to 3.14 Ma.Site U1308 is ∼ 500 m deeper than sites 607 and U1313 and is located to the east of the Mid-Atlantic Ridge, whereas sites 607 and U1313 are in the western Atlantic Basin.The benthic δ 18 O of Site U1308 is similar to sites 607 and U1313 but is of higher resolution (Figs. 3 and 4).Interglacial δ 18 O values are nearly the same for the two sites for the last 1.5 Myr, and slightly lower at Site 607 than U1308 prior to that time.Glacial δ 18 O values tend to be greater at Site U1308 than at sites 607 or U1313.In the latest Pliocene, glacial δ 18 O values show a progressive increase culminating in MIS 100, 98, and 96 (Fig. 4).Following MIS 96, the glacial stages were generally weaker from MIS 94 through 52, with the exception of MIS 82 (2.15 Ma), which was a strong glacial stage in the early Pleistocene (Raymo et al., 1986).Beginning with MIS 52 (1.5 Ma), the amplitude of glacial-interglacial cycles increased because of an increase in glacial δ 18 O and a decrease in interglacial δ 18 O values.From MIS 51 through 25, interglacial δ 18 O was close to the Holocene value of 2.8 ‰ except for MIS 31, 37, and 47, which were particularly strong interglacials.Ruddiman et al., 1986;Raymo et al., 1986).Prominent marine isotope stages are labeled and δ 18 O values corresponding to MIS 5b, 4, and 2 are marked for reference.The orange arrows on the benthic δ 13 C record highlight two trends of decreasing glacial δ 13 C from MIS 22 to 14 and again from MIS 12 to 2 originally described by Raymo et al. (1990).
strong compared to the preceding glacials but not in comparison to MIS 16 and 12, which were the strongest glacial periods recorded for the last 3.2 Myr.The interval from 790 to 480 ka was marked by "luke-warm" interglacials (MIS 19,17,15,13), and interglacial stages became stronger again beginning with MIS 11.
Wavelet analysis of the benthic δ 18 O record indicates a strengthening of 41 kyr power at 2.6 Ma (Fig. 3).In particular, the period between MIS 52 and 36 was marked by an exceptionally pure 41 kyr cycle.Longer periods in the range of 80-100 kyr began to appear at 900 ka and became dominant at 640 ka.

Carbon isotopes
The benthic δ 13 C of Site U1308 is similar to Site 607 although glacial δ 13 C values are generally lower at Site U1308 than Site 607 or U1313 (Figs. 3 and 4), whereas interglacial δ 13 C values are about the same.The lower δ 13 C values during glacial periods may be a consequence of the higher sample density at Site U1308, its slightly greater water depth (∼ 400 m deeper), its location in the eastern basin (Curry and Lohmann, 1983), or the fact that we measured both C. wuellerstorfi and C. kullenbergi.Although C. kullenbergi can possess a shallow infaunal habitat and lower δ 13 C values compared to C. wuellerstorfi (Hodell et al., 2001), paired analyses show no evidence that the δ 13 C of C. kullenbergi is consistently lower than that of C. wuellerstorfi at Site U1308 (Fig. S2).
The first large decrease in benthic δ 13 C occurred in MIS G6 (2.71 Ma), but this event is not observed as strongly at Site 607 or U1313 (Figs. 3 and 4).The lowest δ 13 C values in G6 consist of three measurements of C. wuellerstorfi, which indicates the low values are not anomalous or a consequence of analyzing C. kullenbergi.Persistent decreases in glacial δ 13 C values to < 0 ‰ began with MIS 100 (2.52 Ma).Even during the early Pleistocene, glacial benthic δ 13 C values are often as low as those of MIS 2 (Figs. 3 and 4).Beginning with MIS 52, glacial δ 13 C values are consistently less than values in the last glacial, and are marked by strong 41-cyclicity between MIS 52 and 36.The Site 607 δ 13 C record shows a step-like decrease in δ 13 C beginning with MIS 52 (Raymo et al., 1990).
The lowest δ 13 C values of the Site U1308 record occur in MIS 22 and 12 (Fig. 3).During the last million years, Site U1308 shows the same trends in benthic δ 13 C as Site 607, where minimum glacial δ 13 C values progressively increase from MIS 22 to 14 and again from MIS 12 to 2 (Raymo et al., 1990).The highest δ 13 C values of the entire Site U1308 record occur in MIS 13.
Wavelet analysis of the benthic δ 13 C record shows dominant 41 kyr power in the interval from 2500 to 640 ka followed by a shift to 100 kyr power after 640 ka (Fig. 3).The interval from MIS 52 to 35 has an exceptionally pure 41 kyr cyclicity.

Bulk carbonate δ 18 O
The δ 18 O of bulk carbonate at Site U1308 reflects relative changes in the proportion of biogenic and reworked (detrital or biogenic) carbonate.It is controlled by the input of reworked carbonate but also by changes in the productivity of calcareous microfossils.If biogenic carbonate production is suppressed, then the δ 18 O of bulk carbonate will reflect the isotopic composition of the remaining carbonate, which consists of reworked (older) biogenic carbonate or detrital carbonate, delivered by icebergs.Hodell et al. (2008) showed that lows in bulk δ 18 O are associated with IRD deposition although there are two distinct types of events: one with low carbonate and the other with elevated carbonate concentration.The low-δ 18 O, lowcarbonate events coincide with peaks in density and Si / Sr and are associated with IRD that is rich in silicate minerals.Counts of IRD and foraminifera confirm that the lows in bulk carbonate δ 18 O correspond to lows in foraminiferal counts and peaks in lithic grains, at least for the last several glacial cycles (Fig. S3; Obrochta et al., 2012Obrochta et al., , 2014)).Lowδ 18 O, high-carbonate events are limited to the last 650 kyr and correspond to the lowest bulk δ 18 O values (< −4 ‰) and are associated with Ca / Sr peaks indicative of high concentrations of detrital carbonate in Heinrich layers (Hodell et al., 2008;Hodell and Curtis, 2008).
A potential source of detrital carbonate is IRD originating from the Labrador Sea, which has δ 18 O values averaging −5.6 ‰ ± 1.5 ‰ (Hodell and Curtis, 2008).At Site U1308, the bulk carbonate δ 18 O decreases to the Hudson Strait end member value for Heinrich events 1, 2, 4, and 5 (Hodell et al., 2008).However, the δ 18 O of bulk carbonate during Heinrich events 3 and 6 attains values of only −2 ‰, similar to most of the events prior to 650 ka.Other possible sources of detrital D. A. Hodell and J. E. T. Channell: Mode transitions in Northern Hemisphere glaciation carbonate include limestone and chalk from the British-Irish Ice Sheet (Scourse et al., 2000;Peck et al., 2006) and reworked Cretaceous and Paleogene chalk from northwestern Europe.
The bulk δ 18 O at Site U1308 is marked by long-term, glacial-interglacial changes as well as abrupt events (minima) during glacial periods that coincide with peaks in IRD (Fig. 5).We compare the bulk carbonate δ 18 O with the benthic foraminifer δ 18 O signal because we expect both signals to increase during glacial stages and decrease during interglacial periods.Divergence from this pattern indicates an allochthonous source of carbonate to Site U1308 (Balsam and Williams, 1993;Hodell and Curtis, 2008).We subtract the bulk carbonate δ 18 O from the foraminifer δ 18 O record to emphasize the differences between the two records (Figs. 5  and S4).
Prior to 2.5 Ma, variations in bulk carbonate δ 18 O closely follow the benthic δ 18 O record with higher values during glacials and lower values during interglacials (Fig. S5).Between 2.5 and 1.8 Ma, there are occasional events where bulk δ 18 O is significantly less than foraminifer δ 18 O; for example, during glacial MIS 100 and 98, with a particularly large difference in MIS 82.After ∼ 1.8 Ma, the bulk δ 18 O record is often interrupted by brief minima during glacial periods.The first cluster of low-δ 18 O events occurred during MIS 64, 62, and 60 between ∼ 1.8 and 1.7 ka (Fig. 5).After 1.5 Ma (MIS 50), almost all glacial stages are marked by abrupt decreases in bulk carbonate δ 18 O (Fig. S6).In the interval from 1500 to 650 ka, most of the lows in bulk carbonate δ 18 O are associated with glacial inceptions and/or terminations (Fig. S6).
A pronounced change in the bulk carbonate δ 18 O signal occurs at 650 ka during MIS 16, coincident with the first occurrence of detrital carbonate layers (Heinrich events) derived from Hudson Strait (Hodell et al., 2008).At this time bulk carbonate δ 18 O values reach the detrital carbonate endmember values of −5 ‰ (Fig. 5).After 650 ka, the bulk carbonate δ 18 O minima (some but not all are Heinrich events sensu stricto) occurred relatively late in the glacial cycle and often on glacial terminations (Hodell et al., 2008;Fig. S7).
Wavelet analysis of benthic − bulk δ 18 O shows strong 41 kyr power between 1500 and 650 ka (Fig. 5).Variability in the suborbital (millennial-scale) band also strengthens significantly after ∼ 1.5 Ma.An increase in the precession band (23-19 ka) occurs near 900 ka.In general, the strength of millennial variability appears to be proportional to power in the precessional band.An increase in 100 kyr power occurred at 450 ka.The power spectrum for bulk carbonate δ 18 O contains significant power at ∼ 100, 41, 21-19, 12.8, 9.7, and 8 kyr.(Fig. 6).Some of the higher frequencies may reflect harmonics and/or combination tones of the primary orbital cycles (Hagelberg et al., 1994).

Natural gamma radiation
Natural gamma radiation (NGR) is produced by the radioactive decay of K, Th, and U isotopes, which are contained in clays but can also originate from heavy minerals or lithic grains.Variations in NGR closely follow the benthic δ 18 O record, increasing during glacial stages and decreasing during interglacial periods (Fig. 7).This pattern reflects a greater input of terrigenous sediment during glacial periods and increased carbonate productivity during interglacials.The first significant increase in NGR occurs in MIS 100 (2.52 Ma), consistent with increased delivery of terrigenous sediment to the subpolar North Atlantic at this time (Naafs et al., 2012;Lang et al., 2014).
Wavelet analysis indicates an increase in 41 kyr power beginning at 2.5 Ma and strengthening at 1.5 Ma.Spectra between MIS 52 and 36 are characterized by well-defined 41 kyr cycles, and 41 kyr power dominates until 640 ka, when a quasi-100 kyr cycle emerges in the record (Fig. 7).

Density
In the late Pleistocene, peaks in density are well correlated with lithic grains per gram at Site U1308 (Obrochta et al., 2014).In addition, density varies on glacial-interglacial cycles and increases slowly down-core owing to sediment compaction.The first large density peak occurs just below the hiatus that removed MIS G1 and G2.Strong peaks are recorded in MIS 100 and 98 (Fig. 7), consistent with the widespread delivery of ice-rafted detritus to the subpolar North Atlantic at this time (Shackleton al., 1984;Kleiven et al., 2002;Bailey et al., 2010Bailey et al., , 2012;;Bolton et al., 2010).Density peaks are associated with certain glacial stages including MIS 86,52,46,40,36,34,and 18.The most outstanding feature of the density record are the large peaks beginning at 640 ka in MIS 16 that are associated with Heinrich layers (Hodell et al., 2008;Channell et al., 2012).
Wavelet analysis reveals that density has generally weak power in the orbital band (Fig. 7).Stronger power in the millennial band begins at ∼ 1.8 Ma (MIS 64).

Magnetic susceptibility and grain size
Magnetic susceptibility reflects the total concentration of magnetic minerals, but at this site is dominated by the con- centration of magnetite owing to its high intrinsic susceptibility, whereas κ ARM /κ is a grain size proxy for magnetite.In general, glacial isotopic stages are associated with higher magnetic susceptibility and a tendency towards coarser magnetic grain sizes (i.e., lower values of κ ARM /κ; Fig. 8).Peaks in IRD abundance are similarly marked by abrupt peaks in magnetic susceptibility and coarsening of the magnetic grain size parameter (Channell and Hodell, 2013).IRD transported by icebergs or sea ice from volcanic source areas (e.g., Iceland) is expected to have a disproportionately large effect on magnetic susceptibility.
Magnetic susceptibility begins to increase during glacial periods at 2. 4 Ma with especially high values in MIS 82.Magnetic susceptibility is also high in glacial periods between MIS 74 and 58 between 2 and 1.65 Ma.From MIS 50 (1500 ka) to 16 (650 ka), there is a regular pattern of glacial increases and interglacial decreases in magnetic susceptibility.From 650 ka onward, the Heinrich events are marked by large peaks in magnetic susceptibility.The main feature of the wavelet analysis of magnetic susceptibility is the activation of the millennial band after ∼ 2 Ma.
For κ ARM /κ, MIS 100, 98, and 96 show a muted increase in magnetite grain size consistent with IRD delivery (Fig. 8).MIS 82 is associated with a strong coarsening of magnetic grain size as are MIS 64, 62, 40, and 34.A marked coarsening of magnetic grain size occurs during Heinrich events from 650 ka onwards.Beginning with MIS 19 at 750 ka, there is a distinct decrease in magnetite grain size during interglacials.

Discussion
The discussion is organized around major climate transitions identified at ∼ 2.7, 1.5, 0.9, and 0.65 Ma.The timing of these transitions is not exact and can vary significantly depending upon the proxy considered, presumably because of leads and lags in different components of the ocean-atmosphere system.For example, age estimates for the intensification of Northern Hemisphere glaciation in the late Pliocene range from 3.5 to 2. 4 Ma (Kleiven et al., 2002).A lesser known transition occurred in the middle Pleistocene near 1.6 to 1.5 Ma (Rutherford and D'Hondt, 2000;Hodell and Venz, 2006;Lisiecki, 2014).The onset and end of the MPT have been variably placed between 1.25 and 0.65 Ma (Clark et al., 2006), with an abrupt shift at 0.9 Ma in some records (Elderfield et al., 2012).These transitions represent times of fundamental reorganization of the climate system, and may mark bifurcation points in a dynamic climate system characterized by multiple stable states (Ditlevsen, 2009).

Intensification of Northern Hemisphere glaciation (NHG)
The intensification of Northern Hemisphere glaciation in the latest Pliocene is well documented in North Atlantic sediments.New data from this study are consistent with previous findings at Site U1308 for MIS G8 through 100 (Bailey et al., 2010(Bailey et al., , 2012)).At Site U1308, benthic δ 18 O exceeds 3.5 ‰ for the first time in MIS G6 (∼ 2.72 Ma) and corresponds to a pronounced decrease in benthic δ 13 C (Figs. 3 and 4).
A brief hiatus at Site U1308 removed the record from ∼ 2.6 to 2.65 Ma including MIS G1 and G2 (Fig. 2).The record of natural gamma radiation indicates a significant increase in terrigenous input relative to biogenic carbonate during glacial periods beginning with MIS 100 (2.52 Ma) accompanied by the onset of a distinct 41 kyr cycle (Fig. 7).Increased terrigenous input during glacials beginning at MIS 100 could have been derived from both dust and IRD (Naafs et al., 2012;Lang et al., 2014), although increased magnetite grain size in MIS 100 implies an IRD source (Fig. 8) as the mineral magnetite is not commonly associated with dust.Den-sity peaks occur in MIS 100 (2.52 Ma) and 98 (2.48 Ma) and indicate the occurrence of IRD (Fig. 7).It appears that the widespread occurrence of IRD throughout the North Atlantic did not occur until MIS 100 (Bailey et al., 2010;Bolton et al., 2010;Lang et al., 2014), which underscores the importance of this glaciation relative to prior cold stages.Geochemical provenance studies of IRD carried out for MIS 100 at Site U1308 suggest multiple sources (Bailey et al., 2010(Bailey et al., , 2012)).The early delivery of IRD was mostly from Archean (Greenland) sources, whereas the source shifts to early Proterozoic and Caledonian-age rocks (probably from Scandinavia and North America) during the full glacial conditions of MIS 100 (Bailey et al., 2013).The δ 18 O of bulk carbonate decreases slightly during MIS 100 and 98, but these events are weak in comparison to later glacial periods (Fig. 5).MIS 100, 98, and 96 were strong glacial periods when benthic δ 18 O values at Site U1308 approach those of MIS 4 (Figs. 3 and 4).
Land-based evidence suggests that the Laurentide Ice Sheet advanced to 39 • N at 2.4 ± 0.143 Ma, with the next similarly extensive advance occurring at ∼ 1.3 Ma (Balco and Rovey, 2010).Glacial benthic δ 13 C values at Site U1308 decreased at 2.7 Ma during MIS G6 (Figs. 3 and 4), and the decrease has been interpreted as indicating a decrease in Northern Component Water (NCW; Raymo et al., 1992).More recently, Lang et al. (2016) emphasized the incursion of southernsourced water into the North Atlantic during glacial periods beginning at 2.7 Ma.In the Southern Ocean, the carbon isotope gradient between intermediate and deep water ( 13 C ID ) increased at ∼ 2.75 (MIS G6), marking the development a chemical divide in the Atlantic Ocean between well-ventilated intermediate water and more poorly ventilated deep water (Fig. 9d).Hodell and Venz-Curtis (2006) speculated that this change may signify increased carbon storage in the deep sea and hence a decrease in atmospheric pCO 2 .This prediction appears to be supported by recent paleo-CO 2 reconstructions (Martinez-Boti et al., 2015;Bartoli et al., 2011) and modeling studies (Lunt et al., 2008;Willeit et al., 2015).
The strong glacial triplet of MIS 100, 98, and 96 represented a temporary intensification of glaciation, likely related to a particular set of orbital conditions that included an exceptionally low eccentricity and dampened precession cycle (Maslin et al., 1998).Generally weaker glacials followed from MIS 94 to MIS 52 (1.54 Ma) with the exception of MIS 82 (2.15 Ma), which was a strong glacial stage in the early Pleistocene that coincided with unusually low obliquity.Rohling et al. (2014) identified MIS 82 as the first "deep glaciation" when sea level was lowered below −70 m.Raymo et al. (1986) detected no IRD between MIS 89 and 95 at sites 607 and 609, whereas IRD was detected in MIS 86 and 88 at Site 609 but not at Site 607.
have been similar to MIS 5b, when modeling results suggest the Laurentide Ice Sheet was reduced in size with separate Québec and Keewatin domes (Fig. S8a; Kleman et al., 2013).The Scandinavian and Barents-Kara ice domes were well developed in MIS 5b but did not extend into central Europe or the British Isles (Kleman et al., 2013), and ice volume was comparable between North America and Europe during MIS 5b.Modeling results suggest that the Eurasian Ice Sheet responds dominantly to obliquity forcing (Abe-Ouchi et al., 2013), which may explain the strong 41 kyr cycle in the benthic δ 18 O signal prior to the establishment of very large North American ice sheets.
Benthic δ 13 C values during the early Pleistocene period were often as low as those in the Last Glacial Maximum (LGM), suggesting a shoaling of NCW, increased Southern Component Water (SCW) influence, and/or decreased ventilation (Fig. 3).This is supported by Nd isotope data of Lang et al. (2016), who infer incursions of Southern Component Water during MIS G6, 100, and 96 that were comparable to the LGM.Although glacial ice volume was substantially reduced during the early Pleistocene relative to the LGM, the impact of ice sheets on deep-water circulation was still substantial.The European Ice Sheet may have had a disproportionately large effect on deep-water formation during the early Pleistocene because of its proximity to the Norwegian-Greenland Sea (Fig. 1).Even the relatively small North American ice sheet in the early Pleistocene may have had a strong effect on deep-water circulation because of its proximity to the Labrador Sea.
The magnitude and spatial extent of millennial-scale variability during the late Pliocene and earliest Pleistocene is uncertain.At Site 984 (61 • 25.507 N, 24 • 04.939 W), Bartoli et al. (2006) reported significant millennial-scale variability during glacial stages following the intensification of NHG (2.9-2.8Ma).In contrast, Bolton et al. (2010) found no significant amplification of millennial variability during the glacials of MIS 100, 98, and 96 at Site U1313 (41 • 0.0679 N, 32 • 57.4386 W).They further suggested that the threshold for amplification of millennial variability was not crossed during the late Pliocene and probably not until the MPT.Alternatively, the occurrence of strong millennial variability, as in the late Pleistocene, may have been limited to higher latitudes than Site U1308/U1313 during the late Pliocene-early Pleistocene.
Magnetic susceptibility begins to show an increase in power in the millennial band after 2 Ma (Fig. 8).This change may be related to glaciation in Iceland, where glaciers did not reach sea level until ∼ 2.0 Ma (Einarsson and Albertsson, 1988;Geirsdóttir, 2004).Bulk δ 18 O and density values imply a series of three IRD events associated with MIS 64, 62, and 60 (Fig. 5).This appears to coincide with an increase in magnetic susceptibility at Site 984, south of Iceland, indicating increased delivery of volcanic IRD by icebergs and/or sea ice in glacial periods beginning at 1.8 Ma with MIS 64 (Channell et al., 2002).McIntyre et al. (2001) examined millennial variability at Site 983, south of Iceland, for two periods in the early Pleistocene (1.86-1.93 and 1.75-1.83Ma), including MIS 64 and 70.They found clear evidence for millennial IRD events that recurred approximately every ∼ 2-5 kyr in these two glacial stages.

The 1.5 Ma transition
The climate transition at ∼ 1.5 Ma represented a fundamental change in the mode of glacial-interglacial climate cycles, yet it has received relatively little attention.MIS 52 (∼ 1.54 Ma) marked an important change at Site U1308 as millennialscale variability increased in the mid-latitude North Atlantic.The increased frequency of bulk carbonate δ 18 O events at 1.5 Ma signals the persistent delivery of detrital carbonate to Site U1308 during glacials from MIS 50 onwards (Fig. 5).
We suggest that the lows in the δ 18 O of bulk carbonate between 1.54 and 0.64 Ma were similar to the non-Heinrich IRD layers in the late Pleistocene, which are marked by peaks in lithic grains and low foraminifer abundance.Many of the bulk δ 18 O lows between 1.5 and 0.65 Ma are associated with glacial inceptions and some with terminations (Figs.S5, S6, and S7).These IRD events likely reflect climate-driven changes in the mass balance of ice sheets as a result of advance and retreat of the grounding line at multiple locations in the circum-Atlantic region (Marshall and Koutnik, 2006).
After 1.54 Ma (MIS 52), benthic δ 18 O consistently exceeded 4 ‰ and the interval from 1.54 to 1.2 Ma (MIS 52 to 36) was marked by an exceptionally well defined 41 kyr cycle.We suggest that glacial boundary conditions may have been similar to MIS 4 during this period.During MIS 4, the Laurentide Ice Sheet expanded over Hudson Strait and a high saddle existed connecting the Québec and Keewatin domes (Fig. S7b; Kleman et al., 2013).Full expansion of the Québec Dome, to an extent comparable to the LGM, occurred along the eastern margin of North America during MIS 4 (Kleman et al., 2013), and an ice stream existed in Hudson Strait, thereby supplying detrital carbonate to the North Atlantic during glacial periods, although not by large dynamic surges typical of late Pleistocene Heinrich events.During MIS 4, ice volume was about twice as great in North America as it was in Eurasia.
The size, position, and height of the North American ice sheets have a strong downstream effect over the North Atlantic (Wunsch, 2006;Roberts et al., 2014;Ullman et al., 2014;Zhang et al., 2014), including the position of the winter and summer sea ice limits (Löfverström et al., 2014; Fig. S9).In turn, iceberg drift (and melting) is affected by atmospheric circulation and tends to follow the zero curl of the wind stress.Sea surface temperature (SST) at Site U1313 shows a strong cooling trend beginning at 1.5-1.6Ma (Fig. 10; Lawrence et al., 2010;Naafs et al., 2012).During glacial stages beginning at 1.5 Ma, SST in the subpolar North Atlantic cooled and zonal SST gradients between sites 982 (57.5 • N) and 1313 (43 • N) decreased (Fig. 10).The polar front moved south and winter sea ice extended to the position of Site U1308 during glacial periods.It is unclear whether the increase in glacial IRD after 1.5 Ma represents an increase in production of icebergs or transport and survivability of icebergs out to Site U1308.At Site U1313 to the south, IRD does not begin to greatly increase until ∼ 900 ka (Fig. 10), indicating that changes in IRD abundance are likely diachronous zonally, occurring earlier in the north and later in the south.
On the basis of the benthic δ 13 C of Site 607, Raymo et al. (1990) suggested a significant decrease in the production of North Atlantic Deep Water (NADW) after 1.5 Ma, which is supported by the reanalysis of Lang et al. (2016)  (Fig. 9c).A change in deep-water circulation at 1.6-1.4 Ma is also supported by Nd isotope studies in the North Atlantic (Khélifi and Frank, 2014), which imply that the overflow of deep waters from the Nordic Seas strongly decreased at this time.Bell et al. (2015) suggested that NADW production was strongest in the early Pleistocene between ∼ 2.0 and 1.5 Ma on the basis of high benthic δ 13 C values in the mid-depth southeastern Atlantic (Site 1264).They attributed the change to increased export of dense overflow water from the Nordic Seas into the abyssal eastern Atlantic, which weakened considerably by 1.3 Ma, indicating a reduction in the overflow across the Iceland-Scotland Ridge.
Hodell and Venz-Curtis ( 2006) identified ∼ 1.55 Ma as an important time when the intermediate-deep ( 13 C ID ) gradient increased in the glacial South Atlantic Ocean, indicating increased glacial suppression of NADW and/or reduced ven-tilation of southern-sourced water (Fig. 9d).Lisiecki (2014) also reported changes in Atlantic circulation at 1.5-1.6Ma as indicated by the appearance of glacial δ 13 C gradients between the intermediate and middle-deep Atlantic δ 13 C stacks.Hodell and Venz-Curtis (2006) speculated that the increase in 13 C ID at 1.55 Ma resulted in increased carbon storage in the deep ocean during glacials, and therefore may have been accompanied by lowering of atmospheric CO 2 .On the other hand, boron-isotope-based CO 2 reconstructions do not support a major decrease in CO 2 at 1.5 Ma, although the existing data are low resolution and additional studies are required (Hönisch et al., 2009).
At 1.55 Ma, climate in the polar regions of the North and South Atlantic became synchronized such that both the Arctic and Antarctic polar fronts were marked by glacialinterglacial migrations at a regular pacing of 41 kyr (Hodell and Venz, 1992).Synchronization of the Northern and Southern Hemisphere from ∼ 1.55 Ma onward may indicate that changes in deep-sea carbon storage and CO 2 variations began to play an increasingly important role in glacial-interglacial climate change.This was also a time when dust and iron accumulation began to increase in the subantarctic South Atlantic, suggesting increased iron fertilization and CO 2 drawdown during glacial periods (Fig. 9b; Martinez-Garcia et al., 2011).
MIS 50 also marks a time of increased variability in the millennial band of bulk carbonate δ 18 O and density at Site U1308 (Fig. 5), indicating an increased occurrence of IRD events in the subpolar North Atlantic.Millennial-scale variability was particularly strong between 1.5 and 0.65 Ma at Site U1308.McManus et al. (1999) suggested that the amplitude and frequency of variability in ice rafting and seasurface temperature proxies increases when ice volume is within a critical window defined by benthic δ 18 O values between 3.5 and 4.5 ‰.From 1.5 to 0.65 Ma, the climate system crossed the 3.5 ‰ threshold more often during the 41 kyr world than it did during the 100 kyr world and rarely exceeded the upper threshold of 4.5 ‰ (Fig. 5); thus, the climate system spent more time in the "DO window" (Sima et al., 2004).The benthic δ 18 O threshold may have also been somewhat lower than 3.5 ‰ in the early Pleistocene, if ice sheets flowed more readily than their late Pleistocene counterparts (Raymo et al., 1998;Bailey et al., 2010).
The Site U1308 record suggests that millennial variability was more active in the central North Atlantic from 1.5 Ma, albeit without the glacial dynamics associated with Heinrich events (Hodell et al., 2008).This is supported by results from Site U1385 that demonstrate that millennial variability was a persistent feature on the Iberian margin from 1.5 Ma onwards (Hodell et al., 2015).Birner et al. (2016) showed that the magnitude and pacing of millennial variability during MIS 38 and 40 at Site U1385 was similar to Dansgaard-Oeschger cycles of MIS 3. At Site 983, Raymo et al. (1998) demonstrated millennial-scale variability in MIS 40 and 44 on the basis of proxies of iceberg discharge and deep-water chemistry.
Many of the inferred IRD peaks correspond to low benthic δ 13 C values at Site U1308 during the past 1.5 Myr, suggesting a link between iceberg discharge and weakening of thermohaline circulation (Hodell et al., 2008).Although Site U1308 represents a single site at 50 • N, extension of sea ice to the mid-latitude North Atlantic and changes in Atlantic Meridional Overturning Circulation (AMOC) have been shown to have widespread climate implications (Chiang and Bitz, 2005).

3.11
The 900 ka event Elderfield et al. (2012) deconvolved the benthic δ 18 O record at Site 1123 in the SW Pacific into its temperature and δ 18 O water components by tandem measurement of Mg / Ca and δ 18 O in benthic foraminifera.They inferred a steplike increase in glacial ice volume at ∼ 900 ka.Berger et al. (1993Berger et al. ( , 1994) also found a step-like increase in δ 18 O values beginning at 900 ka during glacials from shallowdwelling planktonic foraminifera on the Ontong Java Plateau in the western Pacific warm pool (Fig. 11).These records suggest that an increase in glacial ice volume across the MPT occurred abruptly at 900 ka during MIS 22, 23, and 24 (Berger, 1993(Berger, , 1994;;Berger and Jansen, 1994;Elderfield et al., 2012).
The MIS 22-24 interval is often considered to be the first 100 kyr cycle because of its similarity with MIS 1-3.The MIS 21-24 interval constitutes a lengthy glacial because MIS 23 was a weak interglacial and the MIS 24-23 transition is not considered to be a full termination (much like the MIS 4-3 transition).Elderfield et al. (2012) proposed that the ice-volume increase during MIS 22-24 may have occurred on Antarctica in response to weak insolation during MIS 23 that suppressed substantial melting of the ice formed in MIS 24.Sea level lowering during MIS 22-24 may have also permitted the advance of marine-based ice sheets onto the continental shelves in the Northern Hemisphere.The onset of widespread glaciation in northern Europe appears to have occurred during MIS 22 in the circum-Baltic region and in Alpine Europe, at the same time as the expansion of widespread lowland glaciation in North America (Head and Gibbard, 2015).Berger and Jansen (1994) suggested that the Svalbard-Barents Sea and Kara Sea ice sheets advanced over the shelf areas after 1 Ma.An over-consolidated section at Site 910 on the Yermak Plateau further supports the grounding of a marine-based ice sheets from Svalbard, and perhaps the Barents Sea, prior to 660 ka (Flower, 1997).
We observe only minor changes at Site U1308 in physical properties or δ 18 O at 0.9 Ma when global ice volume increased.Instead, most physical properties show a transition at 650 ka coincident with the deposition of the first Heinrich layer derived from Hudson Strait (see Sect.Berger et al., 1993Berger et al., , 1994)).Note the abrupt increase in δ 18 O of glacial stages following 900 ka after MIS 22.
ing MIS 22-24, but neither of these events are associated with peaks in detrital carbonate (e.g., increases Ca / Sr) indicating Heinrich layers.
In contrast, Site 982 on the Rockall Plateau (57 • 30.99 N, 15 • 52.00 W; 1145 m water depth) shows three prominent decreases in bulk carbonate δ 18 O during MIS 22-24 (Fig. 12).These δ 18 O lows are associated with decreases in % carbonate and increases in % IRD.The decreases in bulk δ 18 O during MIS 22-24 indicate the delivery of reworked carbonate to Site 982 at 900 ka.At Site 980/981 in Rockall Trough, increases in reworked nannofossil taxa coincide with IRD peaks (Marino et al., 2011).The nannofossils are mostly of Cretaceous age, from the Campanian-Maastrichtian stages, and derived from the Norwegian shelf and/or the northern North Sea-Denmark area (Marino et al., 2011).Reworked specimens of mostly Cretaceous nannofossils also reach a maximum within the 900-700 ka interval in cores from ODP Leg 104 off Norway (Henrich and Bauman, 1994).We suggest that the minima in bulk carbonate during MIS 22-24 at Site 982 were related to expansion of the Eurasian Ice Sheet as sea-level lowering permitted an increase in ice mass on the continental shelves of Scandinavia and the Barents Sea (Berger and Jansen, 1994).The 900 ka event may not have been associated with a major change in the stability of the Laurentide Ice Sheet as there is no evidence for a Heinrich event associated with MIS 22, although the event probably correlates with ice-volume changes in Antarctica and Eurasia, as well as expansion of lowland glaciation in North America.
As noted by Raymo et al. (1997), global benthic δ 13 C records show a pronounced transient decrease at ∼ 0.9 Ma (Fig. 3), representing a perturbation of mean ocean carbon chemistry.Elderfield et al. (2012) attributed this event to changes in deep-water circulation and erosion of organic carbon due to exposure of slope and upper shelf deposits when sea level dropped to −120 m for the first time.Expansion of southern-sourced water into the deep Atlantic is supported by carbon and Nd isotopes that reflect weaker NADW export to the Southern Ocean beginning about 0.9 Ma (Venz and Hodell, 2002;Pena and Goldstein, 2014).Raymo et al. (1997) noted that, over the last 1.5 Ma, the intensity of δ 13 C decreases during glaciations does not directly match the magnitude of δ 18 O increase.For example, there is a progressive increase in glacial benthic δ 13 C values from MIS 22 to 14 that is repeated again between MIS 12 and 2 (Fig. 3), which is not seen in the benthic δ 18 O record.These trends match the magnitude of the terminal IRD peaks found at Site 982 and coincide with transient decreases in benthic δ 13 C at glacial terminations (Venz et al., 1999).The terminal IRD events and δ 13 C decreases during glacial terminations at Site 982 indicate iceberg melting that lowered surface water salinity and thereby reduced Glacial North Atlantic Intermediate Water, resulting in decreased ventilation of the mid-depth North Atlantic.
We suggest that the size and position of the Eurasian ice sheet was a critical factor responsible for the decoupling of benthic δ 13 C and δ 18 O because of the proximity of the ice sheet to source areas of deep-water formation in the Norwegian-Greenland Sea and subpolar North Atlantic  Lisiecki and Raymo, 2005).Note the large millennial decreases in bulk carbonate δ 18 O beginning in MIS 22-24 at 900 ka.
(Fig. 1).Although the Eurasian ice sheet had a relatively small impact on the global δ 18 O of seawater, it may have had a disproportionally large effect on deep-water formation and benthic δ 13 C, causing a decoupling of the response of benthic δ 18 O and δ 13 C signals.
3.12 The 640 ka transition and the emergence of a quasi-periodic 100 ka cycle Undoubtedly the time of greatest change in stable isotopes and physical properties at Site U1308 occurred in MIS 16 at ∼ 650 ka.This change was discussed by Hodell et al. (2008) and marked a pronounced shift in the style and intensity of glacial-interglacial cycles and IRD delivery to the subpolar North Atlantic.Detrital carbonate layers (Heinrich events) first appeared at Site U1308, recorded by large peaks in density and magnetic susceptibility and minima in bulk δ 18 O and κ ARM /κ (Figs. 5 and 8;Hodell et al., 2008;Channell et al., 2012).The results from Site U1308 are supported by results from Site U1313 to the south (Fig. 1), where organic biomarkers indicative of petrogenic compounds derived from Hudson Strait, and an increase in dolomite / calcite ratios, are observed from 0.65 Ma (Fig. 10;Naafs et al., 2013).
As the glacial cycle lengthened, the quasi-periodic cycle of 100 kyr became firmly established at 650 ka (Figs. 3  and 5).Maslin and Ridgewell (2005) have coined the phrase "eccentricity myth" to describe the incorrect attribution of the increase in 80-120 kyr power to orbital eccentricity.Instead, the saw-toothed climate cycles are defined by every four or five precession cycles (Raymo, 1997) or every two or three obliquity cycles (Huybers and Wunsch, 2005;Huybers, 2009), or some combination thereof.Eccentricity modulates the amplitude of the precession cycle and thus may play a role in pacing terminations, but the direct insolation changes resulting from eccentricity are too small to drive the 100 kyr cycle.After 0.9 Ma, the duration of glacial cycles appears to be quantized as multiples of either precession or obliquity (i.e., 80, 100, or 120 kyr).MIS 16 marked the highest benthic δ 18 O values of the entire record, suggesting very cold bottom water temperatures and/or increased ice volume (Fig. 3).MIS 16 and 12 have been referred to as superglacials, indicating continental ice volumes greater than during the LGM.During full glacial periods beginning with MIS 16, it is likely the domes of the North American Ice Sheet coalesced to form a massive unified North American Ice Sheet extending from coast to coast (Bintanja and van der Wal, 2008), similar to conditions during the LGM (Fig. S8c).Hudson Bay would have been covered by a thick ice sheet, and ice streams may have terminated at peripheral ice shelves along the eastern Canadian seaboard (Hulbe et al., 2004).Widespread lowland glaciation was established in northern Europe in MIS 16, together with important expansion of lowland glaciation in North America (Head and Gibbard, 2015).The coalescing of the ice domes in North America permitted the ice sheet to survive subsequent insolation maxima and skip precession or obliquity cycles, thereby lengthening the glacial cycle and transferring power from the precessional and obliquity band to the 80-120 kyr band.
The growth of very large ice sheets also involved fundamental changes in the dynamics of the Laurentide Ice Sheet by introducing instabilities related to processes such as basal melting, isostatic subsidence, and/or drawdown during marine invasion.In addition, the topographic and albedo changes induced by very large ice sheets may have led to a nonlinear response of the climate system (Zhang et al., 2014;Kleman et al., 2013).Large ice sheets charge the freshwater capacitor, permitting the climate system to respond quickly and strongly when the capacitor is discharged.The magnitude of the response was governed by the volume, rate, and location of freshwater addition relative to areas of active deep-water formation (Fig. 1).
IRD events were more frequent, and detrital carbonate layers absent from the record, prior to 650 ka (Fig. 5).This trend is also evident in the NGR record of Site U1304 that shows greater millennial variability prior to 650 ka than afterwards (Fig. 13).As argued previously, millennial variability was more frequent in the 41 kyr world because the climate system spent more time in an intermediate ice-volume state and within the "DO window" (Sima et al., 2004).The decrease in the frequency of IRD events during the past 650 kyr may have also been related to formation of ice shelves, which act as filters of IRD (Alley et al., 2005).Icebergs derived from ice shelves tend to have low debris concentrations because much of the debris is lost by basal melting before the iceberg is calved (Drewry and Cooper, 1981).IRD events after 650 ka occurred less frequently but were greater in magnitude than those prior to MIS 16 (Hodell et al., 2008).Ice-sheet dynamics played a crucial role in the emergence of the saw-tooth pattern of glacial-interglacial cycles in the late Pleistocene.It is not coincidental that Heinrich events (especially terminal Heinrich events) first appeared when the quasi-100 kyr cycle became firmly established as indicated by time series analysis (Mudelsee and Stattegger, 1997;Fig. 3).The development of massive ice sheets during glacial periods beginning at 0.65 Ma (MIS 16) introduced a new type of millennial variability related to episodes of internal dynamical instability of the Laurentide Ice Sheet in the region of Hudson Strait (Hodell et al., 2008).

Co-evolution of millennial and orbital climate variability
The general trend of climate during the Plio-Pleistocene has been one of progressive buildup of larger ice sheets on Northern Hemisphere continents and increased amplitude of glacial-interglacial cycles.Although it is somewhat arbitrary to identify precise change points in the benthic δ 18 O record (Fig. 9), we propose mode transitions in climate evolution during intensification of NHG at approximately 2.7, 1.5, 0.9, and 0.65 Ma.The timing of some of these mode changes is supported by Bayesian change-point analysis of the LR04 stack (Ruggieri, 2012).The benthic − bulk carbonate δ 18 O difference at Site U1308 also shows a trend of increased variability through the Quaternary (Fig. 9a), representing a progressive increase in IRD delivery to the mid-latitudes of the North Atlantic.There is a clear increase in IRD delivery to Site U1308 beginning about 1.5 Ma that coincides with greater amplitude of the benthic δ 18 O signal.At 0.65 Ma, Heinrich events appeared for the first time, introducing a new style of glacial dynamic related to exceptionally large ice sheets and longduration glacial periods (Hodell et al., 2008).
As ice sheets reached the coast, they began to interact with the ocean and affect ocean circulation in areas of deepwater formation such as the Nordic and Labrador seas.As ice sheets continued to grow in size, they stored greater volumes of freshwater on land close to the sources of deep-water formation in the North Atlantic (Fig. 1), increasing the potential of triggering a large climate response.Thus, the magnitude of millennial-scale climate variability may be related to icesheet size and thickness (McManus et al., 1999;Weirauch et al., 2008;Hodell et al., 2008).
We suggest that the variability on orbital and millennial scales may be intrinsically linked and therefore co-evolved during the Quaternary.There is much discussion about the role that millennial variability plays in glacial terminations and inceptions.It is uncertain whether suborbital variability is merely a symptomatic feature of glacial climate or, alternatively, plays a more active role in the inception and/or termination of glacial cycles.For example, strong terminal millennial events may lead to CO 2 degassing in the Southern Ocean through a bipolar seesaw mechanism (Denton et al., 2010;Skinner et al., 2012), thereby hastening deglaciation.Millennial variability may also play a role in glacial ice-sheet buildup by diminishing melting and stabilizing ice sheets during insolation maxima (Timmermann et al., 2010).
Millennial variability acts as an "agitator" of the climate system and can trigger transitions as bifurcation points are approached in a multi-stable dynamical system or act as an activator in an excitable system (Crucifix, 2012).Millennial variability can be considered a form of "noise" on orbital timescales even though it may be derived from deterministic processes.Noise intensity has been found to be important for stochastic resonance in systems containing a subthreshold pacemaker (Perk and Gosak, 2008), which may be important for explaining how relatively weak orbital forcing is amplified by the climate system.The noise level varied on glacial-interglacial timescales and was enhanced during glacial stages and suppressed during interglacials (Hodell et al., 2015), exhibiting characteristics of an excitable system.The interaction of millennial-and orbital-scale climate variability may be an important missing element for explaining the observed patterns of Quaternary climate change.

Caveats
IRD is a tricky proxy to interpret.The mere occurrence of IRD at a particular site does not necessarily imply there was a large climate response to that event.For example, smaller background IRD events are expected to result in a more muted climate response than much larger Heinrich events (Marshall and Koutnick, 2006).The absence of IRD at one site does not necessarily preclude its presence and freshwater forcing elsewhere.Melting of clean icebergs (devoid of IRD) will leave no record of IRD in the sediment yet may result in a significant freshwater forcing.These complications in IRD interpretation have stimulated discussion about the degree to which iceberg discharge is the cause of climate change versus the consequence of stadial conditions (Barker et al., 2015).
Although the benthic δ 18 O threshold concept of millennial variability seems fairly robust (McManus et al., 1999), its exact physical significance remains uncertain because benthic δ 18 O is a function of both ice volume and temperature, and the relative importance of the two has fluctuated during the Quaternary (Elderfield et al., 2012).In addition, the mean δ 18 O of continental ice sheets was likely different (less negative) prior to the late Pleistocene, suggesting the δ 18 O seawater -ice volume relationship has not been constant.The occurrence of millennial variability appears to be related to Northern Hemisphere ice-sheet size as recorded by ben-D.A. Hodell and J. E. T. Channell: Mode transitions in Northern Hemisphere glaciation thic δ 18 O, but the relationship is not entirely straightforward.A critical factor may be how ice growth affects the volume, rate, and location of freshwater discharge to the North Atlantic Ocean relative to the source areas of deepwater formation.In this regard, the growth of the European Ice Sheet may be disproportionately important (relative to its size) for millennial variability because of its proximity to source areas of deep-water formation.
Lastly, our findings are from a single site in the central North Atlantic and may not be representative of the North Atlantic as a whole.Cores from another region may have a different expression of millennial variability in both the magnitude and temporal pattern of change.Indeed, many long records of millennial variability are needed to determine the regional patterns of the time-varying climate signal in the North Atlantic.

Conclusions
Site U1308 provides a 3.2 Myr record documenting the increase in the intensity of Northern Hemisphere glaciation with mode changes at ∼ 2.7, 1.5, 0.9, and 0.65 Ma.The 2.7 Ma transition (MIS G6) marked the appearance of IRD at Site U1308 when benthic δ 18 O first exceeded 3.5 ‰.The event was also associated with a strong decrease in benthic δ 13 C signalling a shoaling of the overturning circulation cell and increased influence of southern-sourced waters in the deep North Atlantic (Lang et al., 2016).Eurasian ice sheets may have had a disproportionally large impact on deep-water circulation during the early Pleistocene because of their proximity to deep-water source areas in the Nordic Seas.The carbon isotope gradient between intermediate and deep water increased at ∼ 2.75 Ma (MIS G6), marking the development a chemical divide in the Atlantic between well-ventilated intermediate water and more poorly ventilated deep water (Hodell and Venz-Curtis, 2006).Increased deep-sea carbon storage resulted in a CO 2 decline (Martinez-Boti et al., 2015) that, together with an exceptional eccentricity minimum, ushered in the glacial-interglacial cycles of the Quaternary.
At Site U1308, the variability in the suborbital (millennialscale) band began to increase at 1.8 Ma and then strengthened significantly after ∼ 1.5 Ma.From 1.5 Ma onward, millennial-scale IRD events became a persistent feature of the Site U1308 record during glacial periods (Fig. 5).The 1.5 Ma transition was also associated with reduction of NCW in the deep North Atlantic during glacials (Raymo et al., 1990), and an increase in vertical carbon isotope gradients between the intermediate and deep ocean, suggesting changes in deep carbon storage (Hodell and Venz-Curtis, 2006).Glacial-interglacial climate changes became synchronized between the subpolar North and South Atlantic at 1.5 Ma (Hodell and Venz, 1992), and CO 2 variations may have assumed an increasingly important role in the synchronization of glacial-interglacial climate change.
No major change is observed at Site U1308 when global ice volume increased at 0.9 Ma; however, at Site 982 (57.5 • N) the MIS 22-24 interval is marked by the first lows in bulk carbonate δ 18 O, indicating delivery of reworked carbonate (Fig. 12), which was likely sourced from the Eurasian ice sheet.We suggest that the 0.9 Ma event may not have involved a major change in the stability of the Laurentide Ice Sheet as there is no evidence for a Heinrich event associated with MIS 22, but rather the ice-volume increase may have been concentrated on Eurasia and Antarctica (Elderfield et al., 2012).The onset of glaciation in the circum-Baltic region and in Alpine Europe is associated with MIS 22 and 20 (Head and Gibbard, 2015).
The time of greatest change in physical properties in the Site U1308 record occurred at 0.65 Ma (MIS 16), when Heinrich events first appeared chronicling a fundamental change in the dynamics of the Laurentide Ice Sheet (Hodell et al., 2008).The intensity of IRD events increased but the frequency decreased at this time, perhaps related to the formation of ice shelves in the Labrador Sea that acted as an IRD filter (Alley et al., 2005).The growth of very large ice sheets on North America and the appearance of Heinrich events in North Atlantic sediments coincided with the emergence of the saw-tooth pattern of δ 18 O change and strong quasiperiodic 100 kyr cycles.
We infer that orbital and millennial variability co-evolved during the Quaternary such that millennial variability generally increased in intensity as ice sheets grew larger in size.Millennial variability provides a source of short-term variability ("agitation") to the climate system that may play an important role in glacial-interglacial climate transitions.The strong link between IRD events and decreases in benthic δ 13 C supports a connection between ice-sheet variability and deep-ocean circulation.Furthermore, the 13 C gradient between intermediate and deep-water increased at the 2.7 and 1.5 Ma transitions (Hodell and Venz-Curtis, 2006), which may reflect an increasingly important role of deep-sea carbon storage and atmospheric CO 2 in glacial-interglacial cycles through the Quaternary.
The Supplement related to this article is available online at doi:10.5194/cp-12-1805-2016-supplement.
Author contributions.Both authors were involved in the planning and execution of IODP Expedition 303, during which IODP Site U1308 was recovered.The authors have contributed equally to data collection, interpretation, and writing of the manuscript.

Figure 1 .
fig01 Figure 1.Location of sites discussed in text (red circles) relative to reconstructed ice-sheet extent in the Northern Hemisphere (continuous solid black line) showing the inferred locations of ice streams (thick black arrows; after Denton and Hughes, 1981) at the Last Glacial Maximum.Yellow stars indicate modern regions of deep water formation and the green star marks the approximate location of intermediate water formation during glacial times.This map emphasizes the point that most of the water stored in the North American and European ice sheets was discharged into the Atlantic Ocean in proximity to areas of deep water formation.Figure modified after Stokes and Clark (2001).

Figure 7 .
Figure 7. Natural gamma radiation (NGR; green) and density (gray) relative to the benthic δ 18 O record of Site U1308.Prominent marine isotope stages are labeled.Top and bottom panels are continuous wavelet transform of the NGR and density records of Site U1308, respectively.Black contours on the wavelet plots indicates the 5 % significance level against red noise.The horizontal dashed white lines represent the orbital periods of precession (∼ 21 kyr), obliquity (∼ 41 kyr), and eccentricity (∼ 100 kyr).

Figure 8 .
Figure 8. κ ARM /κ (proxy for magnetic grain size) and magnetic susceptibility (brown) relative to the benthic δ 18 O record of Site U1308.Prominent marine isotope stages are labeled.Top and bottom panels are continuous wavelet transform of the κ ARM /κ and magnetic susceptibility records of Site U1308, respectively.Black contour on the wavelet plots indicates the 5 % significance level against red noise.The horizontal dashed white lines represent the orbital periods of precession (∼ 21 kyr), obliquity (∼ 41 kyr), and eccentricity (∼ 100 kyr).

Figure 13 .
Figure 13.Natural gamma radiation (NGR) and benthic δ 18 O at IODP Site U1304 on the southernmost Gardar Drift (Expedition 303 Scientists, 2006c; Xuan et al., 2016).Variability in the NGR signal suggests more frequent millennial-scale variability prior to 650 ka (MIS 16) when the climate system spends more time in an intermediate ice-volume state.Magnetostratigraphy is after Xuan et al. (2016).The age model since 1 Ma is based on correlation of the benthic δ 18 O record to LR04, and that prior to 1 Ma on the assumption of constant sedimentation rates between the ages of polarity reversals.