A 1.5-million-year record of orbital and millennial climate variability in the North Atlantic

Climate during the last glacial period was marked by abrupt instability on millennial timescales that included large swings of temperature in and around Greenland (Daansgard–Oeschger events) and smaller, more gradual changes in Antarctica (AIM events). Less is known about the existence and nature of similar variability during older glacial periods, especially during the early Pleistocene when glacial cycles were dominantly occurring at 41 kyr intervals compared to the much longer and deeper glaciations of the more recent period. Here, we report a continuous millennially resolved record of stable isotopes of planktic and benthic foraminifera at IODP Site U1385 (the “Shackleton Site”) from the southwestern Iberian margin for the last 1.5 million years, which includes the Middle Pleistocene Transition (MPT). Our results demonstrate that millennial climate variability (MCV) was a persistent feature of glacial climate, both before and after the MPT. Prior to 1.2 Ma in the early Pleistocene, the amplitude of MCV was modulated by the 41 kyr obliquity cycle and increased when axial tilt dropped below 23.5^∘ and benthic δ¹⁸O exceeded ∼3.8 ‰ (corrected to Uvigerina), indicating a threshold response to orbital forcing. Afterwards, MCV became focused mainly on the transitions into and out of glacial states (i.e. inceptions and terminations) and during times of intermediate ice volume. After 1.2 Ma, obliquity continued to play a role in modulating the amplitude of MCV, especially during times of glacial inceptions, which are always associated with declining obliquity. A non-linear role for obliquity is also indicated by the appearance of multiples (82, 123 kyr) and combination tones (28 kyr) of the 41 kyr cycle. Near the end of the MPT (∼0.65 Ma), obliquity modulation of MCV amplitude wanes as quasi-periodic 100 kyr and precession power increase, coinciding with the growth of oversized ice sheets on North America and the appearance of Heinrich layers in North Atlantic sediments. Whereas the planktic δ¹⁸O of Site U1385 shows a strong resemblance to Greenland temperature and atmospheric methane (i.e. Northern Hemisphere climate), millennial changes in benthic δ¹⁸O closely follow the temperature history of Antarctica for the past 800 kyr. The phasing of millennial planktic and benthic δ¹⁸O variation is similar to that observed for MIS 3 throughout much of the record, which has been suggested to mimic the signature of the bipolar seesaw – i.e. an interhemispheric asymmetry between the timing of cooling in Antarctica and warming in Greenland. The Iberian margin isotopic record suggests that bipolar asymmetry was a robust feature of interhemispheric glacial climate variations for at least the past 1.5 Ma despite changing glacial boundary conditions. A strong correlation exists between millennial increases in planktic δ¹⁸O (cooling) and decreases in benthic δ¹³C, indicating that millennial variations in North Atlantic surface temperature are mirrored by changes in deep-water circulation and remineralization of carbon in the abyssal ocean. We find strong evidence that climate variability on millennial and orbital scales is coupled across different timescales and interacts in both directions, which may be important for linking internal climate dynamics and external astronomical forcing.

1 Introduction

1.1 History of millennial climate variability

Millennial climate variability (MCV) is operationally defined as having a recurrence time between 10³ and 10⁴ years. It excludes variation on orbital timescales but may include harmonics or combination tones of the orbital cycles that have a period of <10 000 years (Berger et al., 2006). MCV is part of the background spectrum of climate variability that follows a power law connecting annual to orbital timescales (Huybers and Curry, 2006). MCV shows closer relationships to Milankovitch cycles than to higher-frequency cycles or oscillations (Huybers and Curry, 2006), and some MCV may result from non-linear coupling of processes operating on orbital timescales (Hagelberg et al., 1994). Because climatic processes are intimately linked across different timescales, documenting the long-term history of MCV is important for understanding its relationship to orbitally forced changes in Quaternary climate.

The first millennial event to be widely recognized in paleoclimate records was the Younger Dryas when a 1300-year-long period of cold climate began at 12 800 yr BP and reversed the general warming trend of the last deglaciation in the Northern Hemisphere (for a review, see Mangerud, 2021). Further study of Greenland ice cores revealed the common occurrence of similar abrupt warming and cooling events during Marine Isotope Stage (MIS) 3 (∼57 to 29 ka). These Dansgaard–Oeschger (D–O) events represent the rapid switching of North Atlantic climate between colder stadial and warmer interstadial states in less than 100 years with a recurrence time of ∼1500 years (Dansgaard et al., 1982). The discovery of such abrupt climate changes in Greenland ice cores in the early 1980s was unexpected because of the great magnitude and rapidity of the temperature change and short recurrence times.

Following the recognition of MCV in Greenland, the search began to see if similar events were recorded in marine sediment cores in the North Atlantic. Marine evidence for D–O events was found in variations in sediment colour and the abundance of the polar foraminifer Neogloboquadrina pachyderma (sinistral) at DSDP Site 609 (Broecker et al., 1990; Bond et al., 1992, 1993). During some of the most extreme stadial events, North Atlantic marine sediment cores were also found to contain layers of ice-rafted detritus (IRD) that are rich in detrital carbonate derived from Paleozoic bedrock underlying Hudson Strait (Heinrich, 1988; Broecker et al., 1992; Hemming, 2004). These so-called “Heinrich events” were attributed to massive discharges of the Laurentide ice sheet to the North Atlantic via Hudson Strait. The D–O cycles are packaged into longer-term cycles (“Bond cycles”) where the amplitude and duration of stadial–interstadial events decrease as climate become progressively cooler until terminating in a Heinrich stadial, which is followed by a large abrupt warming (Bond et al., 1993). The recurrence time of Bond cycles and Heinrich events is on the order of every ∼7–8 kyr, which is longer than D–O events.

MCV, as expressed in Greenland temperature, has a counterpart variation in Antarctic ice cores that is smaller in magnitude and more gradual in nature than the signals found in Greenland. The one-to-one coupling between these events is often explained by changes in interhemispheric heat transport referred to as the thermal bipolar seesaw (Bender et al., 1994; Stocker, 1998; Blunier and Brook, 2001; EPICA Community Members, 2006; WAIS Divide Project Members, 2015). The duration of stadials in Greenland is linearly correlated with the strength of warmings in Antarctica (EPICA Community Members, 2006; WAIS Divide Project Members, 2015). The longer-duration interstadials in Antarctica (Antarctic Isotope Minimum or AIM events) are also marked by rises in atmospheric CO₂ (Ahn and Brook, 2014; Bauska et al., 2021), presumably from decreased stratification and increased overturning in the Southern Ocean (Anderson et al., 2009; Skinner et al., 2010, 2020). On millennial timescales, CO₂ closely tracks Antarctic temperature, with peak CO₂ levels lagging behind peak Antarctic temperature by more than 500 years (Bauska et al., 2021). The magnitude of the CO₂ rise is correlated with the duration of the North Atlantic stadial stage (Buizert and Schmittner, 2015), with a greater CO₂ response during times of prolonged stadial conditions in Greenland, such as those associated with Heinrich events. These longer-lived millennial events represent major reorganizations of the ocean–atmosphere system and have far-reaching effects well beyond the North Atlantic region.

Figure 1 Location of IODP Site U1385 and selected piston (MD95-2042, MD01-2444, MD01-2443) and kasten (MD99-2334K) cores on the Promonotório dos Principes de Avis along the continental slope of the southwestern Iberian margin. The map was made with GeoMapApp (https://www.geomapapp.org/, last access: 15 July 2022) using bathymetry of Zitellini et al. (2009).

A leading hypothesis is that changes in deep-water and/or ocean circulation have played a key role in MCV (for review, see McManus et al., 2004; Alley, 2007; Henry et al., 2016; Lynch-Stieglitz, 2017; Menviel et al., 2020). The Atlantic meridional overturning circulation (AMOC) is sensitive to mode jumps that can be triggered by changes to the surface-water density in North Atlantic source areas of deep-water formation. Climate models of varying complexity have simulated millennial oscillations when forced by freshwater fluxes from melting ice (Stocker and Johnsen, 2003; Ganopolski and Rahmstorf, 2001; Timmermann et al., 2003; Rahmstorf et al., 2005), whereas others have emphasized the role of sea ice (Gildor and Tziperman, 2001; Sevellec and Fedorov, 2015; Li et al., 2005, 2010) and/or ice shelf dynamics (Dokken et al., 2013; Petersen et al., 2013). Some model simulations have shown spontaneous oscillation of the AMOC, even in the absence of deliberate fresh-water forcing (Winton and Sarachik, 1993; Sakai and Peltier, 1999; de Verdiere, 2007; Kleppin et al., 2015). Others have implicated orbitally induced insolation changes or variations in atmospheric CO₂ as (external to the North Atlantic) triggers of MCV (Friedrich et al., 2010; Zhang et al., 2017, 2021; Yin et al., 2021; Vettoretti et al., 2022).

Oxygen isotope records of foraminifera capable of resolving orbital-scale variations are numerous (for a summary of records and resolutions, see Fig. 2 of Ahn et al., 2017), but few long, millennial-resolved records exist to examine the interaction between orbital and millennial components of the climate system. The study of long-term changes in MCV requires long, continuous sedimentary sequences with high sedimentation rates from climatically sensitive areas of the world ocean. Some marine records of MCV exist beyond the last glacial cycle (McManus et al., 1999; Hodell et al., 2008; Oppo et al., 1998; Kawamura et al., 2017; Jouzel et al., 2007; Loulergue et al., 2008; Barker et al., 2011, 2015; Martrat et al., 2007; Margari et al., 2010; Alonso-Garcia et al., 2011; Burns et al., 2019; Gottschalk et al., 2020), but only a few extend beyond 800 ka into the early Pleistocene (Raymo et al., 1998; McIntyre et al., 2001; Birner et al., 2016; Billups and Scheinwald, 2014; Hodell et al., 2015; Hodell and Channell, 2016; Barker et al., 2021, 2022).

Here, we present a 1.5-million-year record of millennial variability in surface- and deep-water properties as recorded by stable isotopes of planktic and benthic foraminifera at IODP Site U1385 (the “Shackleton Site”) located off Portugal in the NE Atlantic Ocean (Fig. 1). The Iberian margin is a well-known location for sediment cores that capture orbital- and millennial-scale variations in North Atlantic climate (Bard et al., 2000; Shackleton et al., 2000, 2004; Martrat et al., 2007; Hodell et al., 2013, 2015). Because of its location in the eastern Atlantic at ∼37^∘ N, the region is sensitive to migrations in the Polar Front but is positioned far enough south that proxies do not saturate under full glacial or interglacial conditions.

The long, millennial-resolved isotope records from Site U1385 provide an opportunity to address several questions about the nature of MCV on orbital and millennial timescales. How common was MCV during older glacial periods of the Pleistocene? Does the nature (intensity, duration, pacing) of MCV change with orbital configuration or climate background state (ice volume, sea level, ice sheet height)? What is the relationship between MCV and longer-term, orbitally driven glacial–interglacial cycles, and how do they interact? How did MCV change across the Middle Pleistocene Transition (MPT) when ice sheets grew larger in size and the amplitude of glacial–interglacial cycles increased? Was the thermal bipolar seesaw mechanism active during older glacial periods of the Pleistocene? What role did millennial variability play in atmospheric CO₂ variations or vice versa?

Figure 2 Comparison of Iberian margin δ¹⁸O records and polar ice cores. Top panel: planktic δ¹⁸O from core MD95-2042 (Shackleton et al., 2000) compared with CH₄ from the West Antarctic Ice Sheet (WAIS) Divide ice core from Antarctica (WAIS Divide Project Members, 2015); Bottom panel: benthic δ¹⁸O in core MD95-2042 compared with the δ¹⁸O record of the WAIS Divide ice core (WAIS Divide Project Members, 2015). Vertical dashed lines are drawn at the abrupt transitions from cold stadials to warmer interstadial conditions in Greenland and are numbered at the bottom of the figure. Note that the phasing of planktic and benthic δ¹⁸O is the same as that inferred from the CH₄ and δ¹⁸O in the WAIS Divide ice core. This pattern has been interpreted as being indicative of a thermal bipolar seesaw.

1.2 The Iberian margin record

The Greenland and Antarctic ice cores provide continuous paleoclimate records to ∼123 000 (NGRIP Project Members, 2004) and 800 000 years (Jouzel et al., 2007) before present, respectively. Beyond the age of the oldest ice, we must rely upon rapidly accumulating marine sediments to document the older history of short-term climate variability in the North Atlantic. Piston cores from the Iberian margin off Portugal contain clear signals of D–O variability in marine sediments (Shackleton et al., 2000, 2004; Pailler and Bard, 2002; Martrat et al., 2007; Margari et al., 2010, 2020). High accumulation rates provide the temporal resolution needed to capture the relatively brief, abrupt temperature changes observed in the Greenland ice core. Shackleton et al. (2000, 2004) demonstrated that each of the D–O events in Greenland is expressed in the Iberian margin planktic δ¹⁸O signal over the last glacial cycle (Fig. 2). In the same sediment core, the benthic δ¹⁸O signal resembles the δD record in Antarctic ice cores (Shackleton et al., 2000, 2004), capturing each of the Antarctic Isotope Maximum (AIM) events (Jouzel et al., 2007). Because the influence of both Greenland and Antarctic millennial events is co-registered in the same sediment core, the phasing can be determined stratigraphically without the usual limitations associated with determining the absolute ages of short-lived climate events. The observed phasing of isotope signals for the last glacial cycle is consistent with the relative changes in temperature between Antarctica and Greenland deduced from the synchronization of ice core records using methane (Fig. 2; Blunier and Brook, 2001; WAIS Divide Project Members, 2015). This pattern has been interpreted as a manifestation of the thermal bipolar seesaw (Stocker and Johnsen, 2003) and can be used to recognize a similar mode of operation of the ocean–climate system in older ice cores (Loulergue et al., 2008) and Iberian margin sediment cores (Margari et al., 2010).

The benthic δ¹³C signal of deep cores from the Iberian margin provides a record of changes in the δ¹³C of deep-water dissolved inorganic carbon (DIC), which varies with changes in deep-water source areas, mixing of water masses, and oxidation of organic matter once the water mass is isolated from the surface ocean. In Iberian margin piston cores, surface cooling is associated with systematic decreases in benthic carbon isotopes, indicating concomitant changes in North Atlantic surface temperature and deep-water circulation (Martrat et al., 2007). Cooling is associated with a shoaling of the Atlantic overturning cell that results in a decreased influence of high-δ¹³C North Atlantic Deep Water (NADW) and an increase of southern-sourced waters with low δ¹³C at abyssal depths in the North Atlantic.

Because of the relative sensitivity of surface and deep-water signals on the Iberian margin to millennial climate change, this area was targeted by the Integrated Ocean Discovery Program (IODP) to extend the record beyond the oldest piston cores from the region. In 2011, five holes were drilled at IODP Site U1385 (the “Shackleton site”) off Portugal, resulting in the recovery of a continuous 166.5 m sequence. A composite section was constructed by correlating elemental data measured by core-scanning XRF at 1 cm resolution in all holes (Hodell et al., 2015). The U1385 record extends to 1.45 Ma (MIS 47) with an average sedimentation rate of 11 cm kyr⁻¹ (Hodell et al., 2013, 2015). The record is mostly complete except for a short hiatus at Termination V that has removed part of late MIS 12 and early MIS 11 (Oliveira et al., 2016).

Figure 3 δ¹⁸O (per mil, Vienna Pee Dee Belemnite (VPDB)) of the planktic foraminifer Globigerina bulloides (a) and mixed benthic foraminifera, mainly Cibicidoides wuellerstorfi, corrected to Uvigerina (b), at IODP Site U1385. Interglacial Marine Isotope Stages (MIS) are labelled with odd-numbered interglacials in (a) and even-numbered glacial stages in (b).

2 Materials and methods

2.1 IODP Site U1385 (“Shackleton site”)

Site U1385 is located very near the position of piston core MD01-2444 (37^∘33.88^′ N, 10^∘8.34^′ W; 2656 m b.s.l. – below sea level; Fig. 1), which consists of a 27 m long sequence representing the last 194 kyr of sediment deposition. Core MD01-2444 can be precisely correlated to Site U1385 on the basis of Ca/Ti measured every 1 cm in both cores (Hodell et al., 2015), thereby providing an equivalent depth (corrected revised metres composite depth – crmcd) in Site U1385 corresponding to each depth in core MD01-2444. Placing MD01-2444 on the Site U1385 depth scale corrects for the well-known effects of stretching and compression that may affect cores recovered with the jumbo Calypso coring system (Skinner and McCave, 2003). Because we did not measure stable isotopes for the upper 23 m of Site U1385 at high resolution, the isotope records presented here consist of a splice between core MD01-2444 (Vautravers and Shackleton, 2006; Margari et al., 2010; Hodell et al., 2013; Tzedakis et al., 2018) and Site U1385 (this study). The Site U1385 record is appended to MD01-2444 at 27.45 m in the piston core, which is equivalent to 26 crmcd in Site U1385, corresponding to an age of ∼194 kyr.

Oxygen and carbon isotope measurements of planktic and benthic foraminifera from Site U1385 were made at an average temporal resolution of ∼200 years for the last 1.45 million years (Fig. 3). The analytical methods were similar to those described by Hodell et al. (2015). For planktic foraminifera, we used the surface-dwelling species Globigerina bulloides from the 250–300 µm size fraction. We interpret the millennial variations in planktic δ¹⁸O of G. bulloides as reflecting variations in sea surface temperature (SST) in the NE Atlantic, which is supported by the strong inverse correlation of planktic δ¹⁸O and alkenone SST data from Iberian margin cores for the past 400 kyr (Martrat et al., 2007). For benthic foraminifera, we used mostly Cibicidoides wuellerstorfi and occasionally other species of Cibicidoides from the >150 µm size fraction. In samples where specimens of Cibicidoides spp. were absent, we used δ¹⁸O of Uvigerina peregrina or Globobulimina affinis. All δ¹⁸O values for each species were corrected to Uvigerina using the offsets suggested by Shackleton et al. (2000) – i.e. +0.64 for Cibicidoides and −0.3 for G. affinis. We recognize that these offsets may vary slightly with time (Hoogakker et al., 2010) but are not large enough to affect the pattern of benthic δ¹⁸O variation.

The water depth of Site U1385 (2578 m b.s.l.) places it under the influence of Northeast Atlantic Deep Water today, but it was influenced by southern-sourced waters during glacial periods. Variations in benthic δ¹⁸O reflect changes in temperature and the δ¹⁸O of deep water bathing the site, which was affected by ice volume on orbital timescales, albeit with such ice volume signals being transported to the core sites on the timescale of ocean mixing (Duplessy et al., 2002; Skinner and Shackleton, 2005; Waelbroeck et al., 2011). Millennial variations in benthic δ¹⁸O are affected by changes in deep-water temperature and by the water mass endmember isotopic compositions (Shackleton et al., 2000; Skinner and Elderfield, 2007; Skinner et al., 2003, 2007). For benthic δ¹³C, we use only the data from the epibenthic C. wuellerstorfi to monitor changes in deep-water ventilation related to changes in deep-ocean circulation and remineralization of organic carbon.

Core-scanning XRF measurements were made every 1 cm in piston core MD01-2444 (Hodell et al., 2013) and in all holes drilled at Site U1385 (Hodell et al., 2015). The Ca/Ti signal was used to correlate among holes and to define a composite spliced section consisting of intervals from Holes A, B, D, and E to form a total length of 166.5 m. The spliced section used in this study consists mostly of Holes D and E, with a few sections taken from Holes A and B to bridge core gaps. All sample depths are given in corrected revised metres composite depth (crmcd) that are corrected for stretching and squeezing caused by coring distortion (Pälike et al., 2005). Theoretically, the same crmcd should be equivalent in all holes, but in practice, the accuracy of the alignment among holes is dependent upon the scale of the correlative features and variability of the Ca/Ti record. We estimate that Ca/Ti features are correlated to the decimetre level or better.

Orbital and millennial variability at Site U1385 is expressed in sediment compositional changes as reflected by elemental ratios (Hodell et al., 2013, 2015). Detrital sediment supply increases relative to biogenic production during cold periods, which is reflected in an increase in Zr–Sr and a decrease in Ca/Ti (Hodell et al., 2015), which are inversely correlated with one another. During the last glacial cycle, increases in Ca/Ti occur during Greenland interstadials, whereas peaks in Zr–Sr mark the stadials, particularly those containing Heinrich events (Channell et al., 2018). Zirconium is mainly derived from zircon, which is a common detrital mineral formed by magmatic and metamorphic processes or the erosion of sedimentary rocks containing zircon. Strontium is highly correlated to Ca, reflecting biogenic carbonate (CaCO₃) because of the incorporation of Sr into biogenic carbonates. We use Sr to normalize Zr, as opposed to Ca, because the signals (counts) are similar and are both measured at the same time (i.e. during the 30 kV scan).

2.2 Chronology

We have updated previous age models of piston core MD01-2444 and IODP Site U1385 (Hodell et al., 2013, 2015) and provide several alternative timescales so users can choose the chronology that is best suited to their specific application. The age models for MD01-2444 include the following (0 to 194 ka): (1) WAIS Divide (WDC2014) by correlation of planktic δ¹⁸O to WAIS methane between 10 and 60 ka; (2) AICC 2012 for MD01-2444 by correlation of benthic δ¹⁸O to δD of EPICA from 60 to 135 ka and using the tie points of Shin et al. (2020) from 135 to 190 ka during MIS 6; (3) a Corchia speleothem chronology is provided for MIS 5 by correlation of planktic δ¹⁸O to the δ¹⁸O of the stalagmite record (Tzedakis et al., 2018).

The age models from MD01-2444 (0 to 194 ka) are combined with those for Site U1385 (>194 ka) to produce the following chronologies: (1) AICC2012 to 800 ka by iteratively correlating millennial events in Site U1385 planktic δ¹⁸O to EPICA CH₄ (gas age) and benthic δ¹⁸O to EPICA δD (ice age); (2) Greenland Synthetic (0–800 ka) by correlation of the planktic δ¹⁸O to Barker et al. (2011); (3) revised LR04 chronology (Lisiecki and Raymo, 2005) based on correlation of Site U1385 benthic δ¹⁸O to the Prob Stack (0 to 1450 ka) (Ahn et al., 2017); and (4) an orbitally tuned timescale by correlation of L^* to the Mediterranean sapropel stratigraphy of the eastern Mediterranean (Konijnendijk et al., 2015). In general, the tuned timescale of Site U1385 compares favourably with LR04 within the estimated error of the chronology, which is ±4 kyr for the past million years and ±6 kyr for the interval from 1.0 to 1.5 Ma (Lisiecki and Raymo, 2005).

The chronology used in this paper is a hybrid model constructed using a combination of age–depth points from MD01-2444 and U1385. The age model is accurate to a precession cycle (∼23 kyr) but cannot provide exact absolute or relative dates for millennial events. This shortcoming limits the reliability of suborbital spectral peaks and estimation of recurrence times of millennial events. Nonetheless, the relative phasing of signals recording different components of the ocean–atmosphere system can be determined stratigraphically without the need for a timescale that is accurate at suborbital resolution. This is particularly important for inferring the phase relationship between planktic and benthic δ¹⁸O, which reflects the interhemispheric leads and lags of the two polar regions.

3 Results

3.1 Defining millennial variability

To identify millennial events, it is necessary to isolate the high-frequency component of the record by eliminating the low-frequency variations related to direct orbital forcing. We experimented with several methods for accomplishing this task, including high-pass filtering, Gaussian smoothing of the record followed by calculation of a residual, and subtracting the planktic and benthic δ¹⁸O values from one another. Although there are subtle differences in detection of millennial events depending on the method and thresholds used, the fundamental identification of millennial events was similar among methods. For simplicity, we settled on a high-pass Butterworth filter of second order with a cutoff frequency starting at $\mathrm{1}/\mathrm{20}$ kyr. The data were interpolated to equal time steps of 0.2 kyr prior to filtering.

Figure 4 (a) The δ¹⁸O record of G. bulloides at Site U1385. (b) High-pass filter of (a) to remove orbital frequencies and to extract suborbital variability. Stadial (blue circles) and interstadial (red circles) events are identified by values that are greater than 1 SD (standard deviation) from the mean. (c) The number of stadial (blue) and interstadial (red) events in non-overlapping windows of 10 000-year duration. (d) Low-pass filter of benthic δ¹⁸O record (black) used to look up δ¹⁸O values for each millennial event. Horizontal dashed black lines correspond to the benthic δ¹⁸O thresholds marking the window of enhanced millennial variability. (e) The number of millennial events is the sum of the stadial and interstadial events in (c). The orange shade indicates times when there are no millennial events per 10 000 years associated with full interglacial stages. Green shade indicates where there are no millennial events per 10 000 years associated with full glacial stages.

We identified stadial and interstadial events using the “findpeaks” function in MatLab by specifying a peak height that must exceed a threshold defined by a multiplier of the standard deviation of the data (e.g. 1σ or 1.5σ) and a minimum peak duration and recurrence time (1 kyr). We varied the parameters so that the algorithm correctly identifies all known D–O events for the last glacial cycle in Core MD01-2444. The same parameters are then applied to identify millennial events for the entire length of the record.

There is some degree of subjectivity involved in identifying millennial events. If the same event is identified in both the planktic δ¹⁸O and Zr–Sr signals (Figs. 4 and 5), we can be confident the event is robust; however, this is not always the case. Not every millennial event in planktic δ¹⁸O has a corresponding change in Zr–Sr, which preferentially records the strongest of the stadial events. Additionally, the planktic δ¹⁸O record can miss stadial events associated with glacial terminations (i.e. terminal stadial event) because the decrease in the δ¹⁸O of seawater from melting ice overwhelms the δ¹⁸O increase expected from cooling. In this case, we rely on the increase in Zr–Sr to recognize the event. Most terminal stadial events are also associated with a minimum in benthic δ¹³C that can be used as an ancillary indicator of these events near glacial terminations. Forthcoming high-resolution measurements of the alkenone SST proxy at Site U1385 will greatly improve the identification of millennial events, especially those associated with terminations (Rodrigues et al., 2017).

Figure 5 Same as Fig. 4 but for Zr–Sr. The blue shade indicates times when the total number of millennial events equals or exceeds three per 10 000 years, which occurs mostly during intermediate glacial states.

We summed the number of millennial events (stadials + interstadials) over a moving non-overlapping window of 10 kyr for both planktic δ¹⁸O and Zr–Sr. Patterns of millennial variability were similar for the two proxies (Figs. 4 and 5). The number of events per 10 kyr interval changes depending upon the choice of start time of the 10 kyr window and whether the analysis is run forwards or backwards, but the fundamental patterns are not substantially altered. The greatest number of millennial events per 10 kyr interval occurred during MIS 3 and during glacial stages of the early Pleistocene from MIS 38 to 46.

3.2 Description of records

Because it is difficult to distinguish millennial events when the Site U1385 record is plotted full scale (Fig. 3), we describe the time evolution of orbital and suborbital variability in the isotope and XRF records for the last 1.45 Ma in ∼200 kyr increments: 0–200 ka (Fig. 6), 200–400 ka (Fig. 7), 400–600 ka (Fig. 8), 600–800 ka (Fig. 9), 800–1000 ka (Fig. 10), 1000–1200 (Fig. 11), 1200–1450 (Fig. 12). We begin with the last 200 kyr because this is the best known period for MCV that can be used as a benchmark for comparison with MCV in the older intervals. Within each interval the record is described oldest to youngest. The records consist of planktic δ¹⁸O, benthic δ¹⁸O, benthic δ¹³C, and Zr–Sr, with stadial events identified by the grey shading. We use a modified version of the isotope nomenclature of Railsback et al. (2015) for marine isotope stages (MIS) of the last million years and the detrital layer stratigraphy of Channell et al. (2012) for Heinrich events.

3.2.1 MIS 1–7a (0–200 ka)

The interval from 0 to 200 ka consists mainly of the record of MD01-2444, which has been described in previous publications (Martrat et al., 2007; Margari et al., 2010, 2014; Hodell et al., 2013). MIS 6 shows a typical pattern of strong MCV at the time of glacial inception following MIS 7a (Fig. 6). Six millennial events are recognized between ∼195 and 155 ka, with a recurrence time ranging from 3 to 7 kyr (Margari et al., 2010, 2014); these also correspond with carbon dioxide maxima (Shin et al., 2020). Minimum benthic δ¹³C values occur at ∼155 ka during event 6vi, which is associated with very cold alkenone SSTs (Margari et al., 2014). MCV becomes more subdued during the full glacial conditions of MIS 6, followed by Heinrich stadial 11 associated with Termination II. MIS 6 shows a clear pattern of decreasing MCV during the glacial cycle, with suppressed variability at the time of peak glaciation. Loulergue et al. (2008), using ice core methane and Antarctic δ¹⁸O, showed a similar pattern of millennial variability, with five interstadial events identified between 190 and 170 ka but only one event identified between 170 and 140 ka. These patterns are also reflected in marine oxygenation reconstructions from the Southern Ocean (Gottschalk et al., 2020). The close similarity in pattern between planktic δ¹⁸O and methane and between benthic and Antarctic ice δ¹⁸O continues throughout the record (Wolff et al., 2022).

Figure 6 0–200 ka: Planktic δ¹⁸O of G. bulloides (black), benthic δ¹⁸O of Cibicidioides spp. (red), obliquity (solid line), precession (dashed line), benthic δ¹³C of C. wuellerstorfi (blue), and Zr–Sr (grey). Odd marine isotope stages (interglacials) are numbered and shaded yellow. Greenland interstadials are numbered, and prominent Heinrich stadials (H) are identified. Even-numbered glacial periods are shaded blue, with stadial events identified by grey vertical bars. Stadials during MIS 6 are numbered after Margari et al. (2010). Terminations are indicated by vertical dashed black lines, and roman numerals have been placed approximately near the mid-point of the deglaciation, although millennial events on the termination often make it difficult to exactly define this point. Horizontal dashed red lines correspond to the benthic δ¹⁸O thresholds marking the window of enhanced millennial variability. Precession and obliquity are plotted such that boreal summer insolation increases in the up direction. The grey shading indicates times when strong MCV is associated with declining or low obliquity, especially associated with glacial inceptions.

Low-amplitude MCV occurs during MIS 5e (Tzedakis et al., 2018) and is followed by three strong stadial events during MIS 5d. MIS 5b is marked by a single prolonged period of stadial conditions. Millennial events DO 20 and 21, documented in the Greenland ice core, are recorded on the transition from MIS 5a to 4. MCV was relatively suppressed during MIS 4, except for a single event (DO 18). The last glacial cycle is unusual in that it is interrupted by a long period of strong millennial variability during MIS 3 followed by a decreased amplitude during the last glacial maximum (MIS 2) between ∼27 and 19 ka (Fig. 2). MIS 2 is terminated, beginning with Heinrich stadial 1, which marks the start of deglaciation. Termination I includes millennial events that occurred during the deglaciation, including the Bølling–Allerød interstadial and Younger Dryas stadial.

3.2.2 MIS 7c–11 (200–400 ka)

The transition from MIS 11 to 10 was marked by strong MCV (Fig. 7), and features in planktic and benthic δ¹⁸O at Site U1385 can be readily correlated to the EPICA ice core records of methane and δD, respectively (Nehrbass-Ahles et al., 2020). Initially, the events are paced ∼5 kyr apart from 400 to 370 ka, and the recurrence time decreases to ∼3 kyr between 365 and 355 ka. MIS 10 culminates in two prolonged Heinrich stadials (10.1 and 10.2) before Termination IV. Low benthic δ¹³C values (<0 ‰) occur at 338 and 352 ka, associated with HS 10.1 and 10.2. MCV is muted during MIS 9a and 9e but is relatively strong in the period between 9a and 9e. MCV resumes during MIS 8, including three Heinrich stadial events (8.1, 8.2, 8.3) prior to Termination III. Minimum δ¹³C values during MIS 8 occur at 272 ka, associated with a very strong cooling event in alkenone SST (Rodrigues et al., 2017). MCV occurs on the transition from MIS 7e to 7d and is relatively suppressed during MIS 7c.

Figure 7 Same as Fig. 6 but for 200–400 ka. Heinrich stadial 8.1, 8.2, 8.3, 10.1, and 10.2 are labelled after Channell et al. (2012). Note strong MCV on MIS 11 to 10 transition culminating in HS 10.2.

3.2.3 MIS 11c–15c (400–600 ka)

Both MIS 15d and 15b contain two strong stadial events, whereas MCV was suppressed during 15a, c, and e (Fig. 8). MIS 14 was a relatively weak glacial by late-Pleistocene standards, and MCV occurred throughout most of the glacial and especially on the MIS 15a–14 transition. MIS 13 shows relatively low variability, with one stadial event in 13c and two near the 13b–a transition. Strong MCV is recorded on the glacial inception of MIS 12, followed by a trend of declining amplitude towards the peak of MIS 12. A minimum in benthic δ¹³C values of < 0 ‰ occurs in the middle of MIS 12 at 455 ka. A short hiatus (∼30 kyr duration) occurs at the transition from MIS 12 to 11, which removed much of Termination V and early MIS 11.

Figure 8 Same as Fig. 6 but for 400–600 ka. A hiatus indicated by dark grey band has removed Termination V and early MIS 11.

Figure 9 Same as Fig. 6 but for 600–800 ka. Note the progressive decrease in the amplitude of millennial events from MIS 17C to the peak of MIS 16 and from MIS 19c to 18.

3.2.4 MIS 15e–20 (600–800 ka)

The end of MIS 20 is marked by a terminal stadial event and a decrease in benthic δ¹³C at 795 ka (Fig. 9). Following MIS 19, strong MCV occurs on the MIS 19–18 transition, with three distinct millennial oscillations paced at ∼5 kyr, which have been interpreted to reflect the second harmonic of precession (Ferretti et al., 2015; Sanchez Goni et al., 2016). MIS 19 is the oldest interglacial recorded in the EPICA Dome C (EDC) ice core, and three consecutive warming events (AIM) occur on the MIS 19–18 transition (Jouzel et al., 2007; Pol et al., 2010); these were also identified in the CH₄ and CO₂ signals (Loulergue et al., 2008; Lüthi et al., 2008). At Site U1385, the phasing between planktic and benthic δ¹⁸O variations during the MCV on the MIS 19–18 transition is similar to that observed during MIS 3, suggesting an active bipolar seesaw (Fig. 2). The phasing between methane and δD in the EPICA ice core is difficult to determine because of large uncertainties in gas-age–ice-age offsets and possible diffusion in the deepest part of the ice core.

MIS 18 consists of two distinct glaciations separated by a long interstadial period that is punctuated by a stadial event in the middle at 730 ka. Millennial variability decreases throughout the glacial inception towards the first glacial peak, associated with a decrease in benthic δ¹³C at 742 ka. The second glacial peak of MIS 18 is marked by a very strong decrease in benthic carbon isotope values at 717 ka, associated with Terminations VIIIA.

MIS 16 shows a trend of decreasing amplitude of MCV through the glacial cycle where the variability is greatest on the MIS 17–16 glacial transition and diminishes towards the peak glacial conditions of MIS 16. Strong stadial events associated with Heinrich events 16.1 and 16.2 are suspiciously absent near Termination VII, perhaps indicating the presence of a previously undetected hiatus.

3.2.5 MIS 21–27 (800–1000 ka)

The pattern of increased MCV associated with the transitions from interglacial to glacial stages continues with MIS 27–26 (Fig. 10). MIS 26 and 28 were relatively weak glacials and were marked by strong millennial variability. The interval from MIS 25–21 is often compared with MIS 5–1 because of the similarity of weak interglacial MIS 23 to MIS 3. MIS 25–21 is sometimes erroneously described as the first “100 kyr cycle”, but it consists of two obliquity cycles (Bajo et al., 2020). The MIS 24–23 transition (TXI) was an incomplete (skipped) deglaciation, thereby lengthening the duration of glacial conditions to ∼80 ky. The pace of millennial events is faster for the MIS 25–24 transition than for some other glacial inceptions. Strong MCV is evident throughout MIS 24, and unlike MIS 3, MCV is relatively suppressed during MIS 23, which contains a single strong millennial event at 919 ka. This pattern is different from the last glacial cycle, when MCV was suppressed during MIS 4 and enhanced during a significant portion of MIS 3.

Figure 10 Same as Fig. 6 but for 600–1000 ka, which includes the “900 kyr” interval that clearly consists of two obliquity cycles with strong MCV on the MIS 25–24 transition.

Glacial ice volume increased (Elderfield et al., 2012) during MIS 25–21, and major changes occurred in deep-ocean circulation (Pena and Goldstein, 2014) and carbon cycling (Thomas et al., 2022). Minimum benthic δ¹³C values occurred at 878, 898, and 925 ka, which are among some of the lowest values found in the deep North Atlantic during the Quaternary (Raymo et al., 1997; Hodell and Channell, 2016). MIS 21 has multiple substages and consists of four warm periods that are spaced about 10 kyr apart; these have been interpreted as the second harmonic of precession (Ferretti et al., 2015).

3.2.6 MIS 28–36 (1000–1200 ka)

Although the timing of the MPT is ambiguous and dependent upon the proxy signal and method used to define the shift in the dominant period from 41 kyr to longer periods, the lengths of some of the glacial–interglacial cycles appear to increase from ∼1200 kyr as the shape of the benthic δ¹⁸O signal becomes less symmetrical and takes on a more sawtooth waveform (Broecker and van Donk, 1970), reflecting a slower rate of ice sheet growth and a faster rate of decay. For example, MIS 35–34–33 has a different duration and shape than previous glacial cycles – it is drawn out with ∼90 kyr between TXVI and TXIV. MIS 35 was an exceptionally long interglacial (Shackleton et al., 1990). Very strong MCV occurred on the MIS 35–34 transition (Fig. 11), consisting of four prominent events that are paced about 5–6 kyr apart. The abrupt warming events that mark the start of interstadials coincide with minima in benthic δ¹⁸O, indicating that the phase relationship is similar to that observed during MIS 3 between Greenland and Antarctica (Fig. 2), which is a pattern indicative of the bipolar seesaw. The benthic δ¹³C mirrors the planktic δ¹⁸O record, with strong decreases in benthic δ¹³C associated with each of the stadial events.

Figure 11 Same as Fig. 6 but for 1000–1200 ka. Note the very strong MCV on the MIS 35–34 transition associated with declining obliquity.

The stadial events become progressively colder during MIS 34, culminating in the terminal stadial event that occurs near the MIS 34–33 transition (TXV). This stadial is marked by a strong decrease in benthic δ¹³C. Benthic δ¹⁸O begins to decrease first, while planktic δ¹⁸O remains high (cold) and benthic δ¹³C low. This is the same phasing as that observed during Termination I on the Iberian margin, when the Southern Hemisphere begins to warm at ∼18 ka while the North Atlantic remains cold and NADW shoals (Skinner and Shackleton, 2006).

Millennial events occur within MIS 33 and on the glacial inception of MIS 32, with a strong terminal stadial event associated with MIS 32–31 (TXIV). MIS 31 (1094–1062 ka) was also a relatively long and strong interglacial (Oliveira et al., 2017). MCV occurs on the 33–32, 31–30, and 29–28 glacial inceptions, and each is associated with declining obliquity.

3.2.7 MIS 37–47 (1200–1440 ka)

The period from ∼1250 to 1550 ka (MIS 52) in the early Pleistocene was a time when glacial and interglacial cycles, as recorded by benthic δ¹⁸O, were dominated by a 41 kyr period corresponding to variations in Earth's obliquity, although precession was still significant (Liautaud et al., 2020). Similar to the last glaciation and Holocene, MCV is enhanced during glacial periods and suppressed during interglacial stages, exhibiting a threshold response (Fig. 12). MCV increases when obliquity drops below a threshold value of 23.5^∘, which corresponds to a benthic δ¹⁸O threshold of ∼3.8 ‰ (corrected to Uvigerina). Importantly, and unlike late Pleistocene glaciations after the MPT, MCV persists throughout most of the glacial part of the cycle. Many of the increases in planktic δ¹⁸O (stadial events) are associated with coeval decreases in benthic δ¹³C, indicating a link between North Atlantic surface climate and deep-water circulation.

Figure 12 Same as Fig. 6 but for 1200–1440 ka, when glacial–interglacial cycles were occurring every 41 kyr. Note the increase in MCV associated with obliquity lows.

4 Discussion

4.1 Millennial variability in planktic δ¹⁸O

Because of the great similarity between Greenland δ¹⁸O and Iberian margin planktic δ¹⁸O signals for the last glacial cycle, we interpret this proxy of surface temperature as an indicator of MCV in the North Atlantic. XRF records of Site U1385 provided the first evidence that MCV was a persistent feature of the climate system for at least the past 1.45 Ma (Hodell et al., 2015), which is confirmed by comparison of the planktic δ¹⁸O and Zr–Sr signals (Figs. 4 and 5). The first-order pattern is that MCV is enhanced during glacial stages and diminished during each of the full interglacial stages (see shading in Fig. 4), which is consistent with the relative stability of Holocene climate relative to the last glacial period in the Greenland ice core and with other paleoclimatic results (McManus et al., 1999; Barker et al., 2021; Sun et al., 2021; Kawamura et al., 2017). A repeated pattern is that the end of each interglacial stage is marked by the onset of strong MCV that continues through the period of glacial inception when ice sheets expand on North America and Eurasia. MCV associated with glacial inceptions have generally longer recurrence times than D–O events, varying between 3 and 8 kyr. A few glacial cycles of the late Pleistocene show a clear pattern of decreasing amplitude of MCV from the glacial inception towards peak glacial conditions (e.g. MIS 6, 10, 12, and 16; Figs. 6–9), giving rise to a sawtooth shape. The pattern of MCV evolves from longer, stronger interstadials to shorter, weaker interstadials as climate becomes progressively cooler during the glacial cycle. In fact, it is MCV that is partly responsible for the unevenly spaced teeth in the sawtooth pattern of interglacial-to-glacial transitions. The last glacial cycle is unusual in that the low MCV during MIS 2 and 4 is interrupted by a period of strong variability during MIS 3. Such D–O-type MCV has a short recurrence time (1.5–2 kyr), which is also found during early Pleistocene glaciations prior to 1.25 Ma (Birner et al., 2016). Almost all glacial periods end with a strong terminal stadial event that marks the start of deglaciation, with some terminations containing additional millennial events during deglaciation (e.g. Bolling–Allerod and Younger Dryas oscillations).

4.2 Millennial variability in benthic δ¹⁸O

Unlike planktic δ¹⁸O, Shackleton et al. (2000, 2004) demonstrated that variations in the benthic δ¹⁸O signal of piston cores from the Iberian margin closely follow the δD record of Antarctic ice cores for the last glacial period (Fig. 2). The Site U1385 record indicates that this similarity extends for at least the last 800 kyr (Fig. 13; Nehrbass-Ahles et al., 2020). The reason for the similarity of Iberian margin benthic δ¹⁸O and Antarctic temperature is not entirely clear (Skinner et al., 2007). Shackleton et al. (2000) originally proposed that the millennial oscillations in benthic δ¹⁸O during MIS 3 reflected changes in the δ¹⁸O of seawater, caused by ice volume variations of the order of 20–30 m of sea level equivalence (Siddall et al., 2008). An alternative explanation is that millennial variations in benthic δ¹⁸O reflect temperature changes of deep water. In this case, the large variations in Antarctic air temperatures are damped by the thermal mass of the deep ocean and translate into small changes in benthic δ¹⁸O, reflecting temperature changes in the source areas of deep-water formation around Antarctica. The similarity of deep-water temperature estimated by Mg–Ca at ODP Site 1123 in the South Pacific and of Antarctic temperature (Elderfield et al., 2012) supports this interpretation, as does the emerging but sparse evidence for similarity between mean ocean temperature and Antarctic temperature (Haeberli et al., 2021). Surface temperatures in the high-latitude Southern Ocean may have been important for regulating deep-ocean heat content, which has implications for deep-ocean circulation and CO₂ storage (Jansen, 2016). Skinner et al. (2007) measured benthic Mg–Ca and δ¹⁸O in core MD01-2444 during MIS 3 and concluded that the benthic δ¹⁸O record cannot be interpreted as a unique proxy of either deep-water temperature or ice volume and must contain a significant local hydrographic component related to the mixing of endmember water masses from the North Atlantic and Southern Ocean, which have different δ¹⁸O values. This is further supported by similar results from the deep Southern Ocean, where benthic δ¹⁸O exhibits a similar (but not identical) pattern to that observed on the Iberian margin (Gottschalk et al., 2020), and deep-water temperatures again appear to have decreased during HS4, consistent with enhanced convection contributing to Antarctic warmth and atmospheric CO₂ rise (Skinner et al., 2020; Menviel et al., 2014). In all cases, a global glacioeustatic signal would only be transported around the globe on a timescale that is consistent with ocean transport and mixing (i.e. centuries to millennia) (Primeau and Deleersnijder, 2009), which would oppose any proposal of benthic δ¹⁸O tracking global ice volume in synchrony (Gebbie, 2012). Indeed, this is demonstrated by the phasing of benthic δ¹⁸O, Antarctic temperature, mean ocean temperature, and sea level on the last deglaciation (Skinner and Shackleton, 2005; Baggenstos et al., 2019).

Figure 13 Comparison of benthic δ¹⁸O from MD01-2444 and Site U1385 on the Iberian margin (red, this study) and δD from EPICA Dome C ice core, Antarctica, for the last 800 kyr (black, Jouzel et al., 2007). Timescale is AICC2012 (Veres et al., 2013). Interglacial marine isotope stages are numbered and shaded yellow.

As in the latest Pleistocene, stadial events are associated with decreases in benthic δ¹³C for the past 1.45 Ma, suggesting that surface coolings in the North Atlantic were associated with perturbations of deep-water ventilation and carbon storage in the deep Atlantic (Martrat et al., 2007; Shackleton et al., 2000; Skinner et al., 2007). Low δ¹³C values are associated with each of the glacial terminations, when δ¹⁸O decreases, and in some cases, the low δ¹³C values are prolonged and extend into the early part of the interglacial period (Hodell et al., 2009; Galaasen et al., 2014, 2020).

4.3 Phasing of planktic and benthic δ¹⁸O and the bipolar seesaw

Because of the similarity of the planktic and benthic oxygen isotope records to Greenland and Antarctica, respectively, Shackleton et al. (2000) suggested that the relative phasing of interhemispheric climate change could be assessed in the depth domain using a single core from the Iberian margin. The lead of millennial-scale benthic over planktic δ¹⁸O in piston cores from MIS 3 (Shackleton et al., 2000; Skinner et al., 2007) is observed throughout the Site U1385 record (Fig. 14). This apparent lead of the benthic over the planktic δ¹⁸O is more likely the consequence of the different shapes of the benthic (rectangular) and planktic (triangular) signals (Hinnov et al., 2002). Rather than a direct lead–lag relationship between the polar regions, the thermodynamic bipolar seesaw model predicts an anti-correlation between Greenland temperature and the rate of change of Antarctic temperature, with the abrupt warmings in Greenland leading the Antarctic cooling onset by about 200 years (WAIS Divide Ice Core members, 2015). The consistent phase relationships between planktic and benthic δ¹⁸O during some millennial events for the past 1.45 Ma suggest that the oceanic bipolar seesaw was a robust feature of the interhemispheric climate system despite differences in climate background state. For example, the phasing is the same during glacial inceptions as it is during deglaciations and intermediate ice volume states such as MIS 3. Millennial variation in AMOC and the thermal bipolar seesaw represent mechanisms by which MCV can be propagated from the North Atlantic to the broader climate system.

Figure 14 Examples of the phasing of millennial benthic and planktic δ¹⁸O variability in depth domain: (a) MIS 1–5; (b) MIS 9–11; (c) MIS 23–25; (d) MIS 33–35; and (e) MIS 37–39. The vertical dashed lines mark the rapid warmings (decreases) in the planktic δ¹⁸O record (black) and decreases in benthic δ¹⁸O (red). The decrease in benthic δ¹⁸O occurs prior to the decrease in planktic δ¹⁸O, which is similar to the phasing observed during MIS 3 (a and Fig. 2).

4.4 Orbital modulation of MCV

To test for amplitude modulation of millennial variability by orbital cycles, we follow the approach of Hinnov et al. (2002), who analysed the MD95-2042 record for the last 100 kyr. We examine the power spectrum of the planktic δ¹⁸O and Zr–Sr records after applying a Taner filter and Hilbert transform. Bandpass filtering was performed on evenly resampled (0.2 kyr) time series using a Taner filter centred on 0.55±0.45 with a roll-off rate = 1×10¹², which has better leakage suppression outside the stop band compared to the Butterworth filter. The instantaneous amplitude of the modulating signal was calculated by Hilbert transformation. The presence of significant orbital frequencies in the power spectrum of the Hilbert transform indicates orbital modulation of the amplitude of MCV, and the evolutive spectra show how the orbital modulation of MCV has changed through time (Figs. 15 and 16).

Figure 15 Evolutive power spectrum of the amplitude modulation of planktic δ¹⁸O as estimated from a Taner filter (TF) centred at 0.55±0.45 and Hilbert transformation (HT) of the time series. Sliding window of ∼300 kyr with time domain zero padding and a step equal to the sampling rate of the time series (∼0.2 kyr). Evolutionary spectra created with Acycle software in MatLab (Li et al., 2019).

Figure 16 Same as Fig. 15 but for Zr–Sr.

The power spectra of planktic δ¹⁸O and Zr–Sr are similar and support orbital modulation of the amplitude of the millennial band by Earth's orbital parameters (e.g. 19, 23, 28, 41 kyr). The 41 kyr obliquity dominates the modulation of MCV between 1450 and ∼900 ka, with a weak precession component (Figs. 15 and 16). At ∼900 ka, power develops at ∼28 kyr, and precession strengthens, especially in the Zr–Sr record (Fig. 16). The 28 kyr cycle is not entirely unexpected because it has been widely reported in late Pleistocene ice core and marine δ¹⁸O records (Huybers and Wunsch, 2004; Lisiecki and Raymo, 2005; Lourens et al., 2010). We note that the theoretical obliquity signal contains secondary peaks at ∼29 kyr and 54 kyr (Laskar et al., 2004), but their spectral power seem too weak to be of any direct climatic significance. Instead, the 28 kyr cycle has been interpreted by Lourens et al. (2010) as resulting from the sum frequencies between the 41 kyr cycle and its multiples of 82 kyr (i.e. $\mathrm{1}/\mathrm{82}+\mathrm{1}/\mathrm{41}=\mathrm{1}/\mathrm{27.3}$) and 123 kyr (i.e. $\mathrm{1}/\mathrm{123}+\mathrm{1}/\mathrm{41}=\mathrm{1}/\mathrm{30.8}$). However, the 28 kyr power could also result from the difference of frequencies between multiples of the 41 kyr cycle and the main precession components (e.g. $\mathrm{1}/\mathrm{21}-\mathrm{1}/\mathrm{82}=\mathrm{1}/\mathrm{28.2}$). Liebrand and de Bakker (2019) applied bispectral analysis techniques to the LR04 benthic δ¹⁸O signal and showed that a large part of the precession (spectral) energy could have been transferred to the lower frequencies of obliquity and its multiples in the course of the Quaternary, and especially during the MPT. In this respect, the presence of a strong precession signal in both the planktonic δ¹⁸O and Zr–Sr records of U1385 could be partially responsible for the occurrence of the ∼28 kyr beat, but additional bispectral analyses is required to further unravel these energy transfers. At ∼650 ka, the 41 and 28 kyr power of obliquity declines substantially, and the spectrum is marked by lower-frequency power (80–120 kyr), which is difficult to interpret in the evolutive spectrum because of the relatively short window size (∼300 kyr). This may reflect an increase in eccentricity modulation of MCV or modulation by multiples of the obliquity and precession cycles or a change in the non-linear energy transfer between orbital components across the MPT (Liebrand and de Bakker, 2019).

The early Pleistocene interval from 1200–1440 ka provides strong evidence for a relationship between the occurrence of MCV and obliquity (Fig. 12). MCV increases during times of low obliquity, displaying a threshold response such that it increases each time the obliquity drops below ∼23.5^∘. The relationship between MCV and obliquity is not a result of orbital tuning, as the chronology was derived by tuning the colour (L^*) to precession and is independent of obliquity (Hodell et al., 2015). It is uncertain, however, if the increased millennial variability is the direct (i.e. fast-acting) result of obliquity through its effect on the mean insolation (Berger et al., 2010) and mainly on the total summer insolation at high latitudes (e.g. lowered insolation, colder temperature, sea ice expansion) or on the meridional insolation gradient. Alternatively, low obliquity secondarily leads to increased ice volume, raised ice sheet height, and lowered sea level (see Discussion, Sect. 4.5) – all of which have been proposed as a threshold for MCV (McManus et al., 1999; Zhang et al., 2014).

Some modelling experiments have demonstrated increased MCV during times of low obliquity in the absence of freshwater forcing (Friedrich et al., 2010; Brown and Galbraith, 2016; Galbraith and de Lavergne, 2018). The obliquity threshold observed for the early Pleistocene is highly suggestive of a non-linear system that is influenced by orbital cycles. For example, sea ice expansion during times of low obliquity may provide strong albedo feedback amplification, resulting in a non-linear response (Tuenter et al., 2005). As the mean position of the sea ice edge expands to lower latitudes, the region of deep-water formation moves from the Norwegian–Greenland Sea to south of Iceland, shifting the AMOC with respect to the mean atmospheric precipitation field where precipitation exceeds evaporation, thereby making the system less stable (Sevellec and Fedorov, 2015; Friedrich et al., 2010).

The relationship between low obliquity and enhanced MCV persists after 1.2 Ma and is expressed as increased millennial variability associated with the transitions from interglacial to glacial stages, which are always associated with declining obliquity (Tzedakis et al., 2012; O'Neill and Broccoli, 2021). In view of this, Tzedakis et al. (2012) proposed that the end of interglacials could be defined as 3000 years (kyr) before the reactivation of MCV at the time of glacial inception. Low obliquity is important for controlling ice accumulation at the start of a glaciation because ice growth begins at high latitudes (and altitudes) where the effect of obliquity on summer insolation is strongest (Vettoretti and Peltier, 2004). Lower obliquity decreases the summer insolation at high latitudes, reduces seasonality, and strengthens the insolation gradient between low and high latitudes, thereby increasing the meridional heat and moisture flux to the high latitudes (Mantsis et al., 2011). The increased heat transport does not balance the direct cooling effects of obliquity through reduced insolation at high latitudes. Increased moisture transport towards the poles provides the fuel needed for growing ice sheets (Vimeux et al., 1999; Raymo and Nisancioglu, 2003). Reduced temperature and increased moisture flux are the two ingredients needed for rapid ice sheet growth during glacial inceptions. At glacial inceptions, precession and atmospheric CO₂ play secondary roles that may reinforce or delay increased ice accumulation, depending on CO₂ concentration and the phasing of precession and obliquity (Vettoretti and Peltier, 2004).

Modelling studies suggest that orbital forcing may play a more direct role in the onset of MCV at the end of interglacial periods. Using LOVECLIM1.3, Yin et al. (2021) found a threshold response to decreasing summer insolation related to both precession and obliquity. When summer insolation falls below a critical value, a strong, abrupt weakening of AMOC is triggered as sea ice expands in the Nordic and Labrador seas. The transition into a cooler mean climate state is accompanied by high-amplitude temperature variations lasting for several thousand years (Yin et al., 2021).

Zhang et al. (2021) used a fully coupled climate model and found that changes in Earth's orbital geometry can directly affect MCV during intermediate glacial states (e.g. MIS 3). Both obliquity and precession play a role in AMOC stability (Zhang et al., 2021; Yin et al., 2021) in terms of its effect on mean insolation at high latitudes and in eccentricity-modulated precession in terms of its low-latitude effect on the subtropical hydrologic budget and salinity of the North Atlantic basin.

MCV can also result from orbital forcing that is expressed as subharmonics and combination tones of the primary orbital cycles. Using bispectral analysis, Hagelberg et al. (1994) demonstrated that approximately a third of the variability in the frequency band ranging from $\mathrm{1}/\mathrm{15}$ to $\mathrm{1}/\mathrm{2}$ kyr originates from the transfer of spectral energy from the lower-frequency Milankovitch band (see also Liebrand and de Bakker, 2019). A case in point is the 11 and 5.5 kyr cycles found in MIS 21 and 19, respectively, that have been attributed to the second and fourth harmonics of the primary precession cycles (Ferretti et al., 2015; Sanchez Goni et al., 2016). Berger et al. (2006) suggested that the double maximum that occurs in daily irradiation at tropical latitudes includes a suborbital insolation forcing at 11 and 5.5 kyr periods related to precession harmonics.

4.5 State dependence of MCV

Orbital changes may influence MCV directly through fast processes (e.g. sea ice) or more indirectly through slow changes in ice sheet configuration (volume or height) and sea level. On the basis of a 500 kyr long record of ice-rafted detritus and summer SST from Site 980 at 55^∘ N in the Rockall Trough, McManus et al. (1999) suggested that MCV was enhanced during times of intermediate ice volume, as defined by a window or “sweet spot” when MCV was most active during times of intermediate glacial states (Sima et al., 2004; Galbraith and de Lavergne, 2018). MCV is suppressed under full interglacial conditions and during some peak glaciations. The concept of increased MCV during times of intermediate ice volume is supported by observations from the last glacial cycle, when MCV was relatively suppressed during MIS 2 and 4 and strong during MIS 3. MCV was also frequent during glacial periods of the early Pleistocene between 1.45 and 1.25 Ma, when glacial benthic δ¹⁸O values fell entirely within the millennial window (Fig. 5). After 1.25 Ma, the benthic δ¹⁸O threshold is crossed slowly during glacial inceptions and more quickly at glacial terminations (Sima et al., 2004) with some, but not all, full glacial periods marked by reduced MCV.

We tested whether there is a statistically significant tendency for millennial events to occur within a certain range of benthic δ¹⁸O values at Site U1385. The findpeaks algorithm in MatLab returns the age of each event identified, which is then used to look up its corresponding benthic δ¹⁸O value. The δ¹⁸O values are concatenated to form a subpopulation of benthic δ¹⁸O values corresponding to millennial events; this is compared with the full population of benthic δ¹⁸O values (Fig. 17a and b). A two-sample Kolmogorov–Smirnov (K–S) test is used to evaluate if the two populations are from the same or different continuous distributions and whether the tail of the millennial subpopulation distribution is smaller than the full population of benthic δ¹⁸O values.

Figure 17 Probability density estimate of benthic δ¹⁸O for all values (black) and for interstadial (red) and stadial (blue) events for planktic δ¹⁸O (a) and Zr–Sr (b). Vertical dashed lines represent benthic δ¹⁸O threshold values for MCV that define the millennial window.

For millennial events identified in both planktic δ¹⁸O and Zr–Sr, the millennial benthic δ¹⁸O population was significantly different from the full population at 95 % confidence, and the tail of the millennial populations was significantly smaller than that of the full δ¹⁸O population. Millennial events are clearly less frequent at the low (warm) end of the benthic δ¹⁸O distribution, suggesting reduced MCV during full interglacial periods. We estimate the lower benthic δ¹⁸O threshold to be ∼3.8 ‰ for both planktic δ¹⁸O and Zr–Sr (note that 0.64 ‰ must be subtracted from this value to convert to the Cibicidoides scale; Fig. 10c and d). The δ¹⁸O threshold for MCV may differ depending on the record and proxy used to identify millennial variability (IRD, SST, planktic δ¹⁸O, etc.) and may be non-stationary through time. For example, Bailey et al. (2010) suggested that the δ¹⁸O threshold for the late Pliocene (MIS G4 at ∼2640 ka and MIS 100 at ∼2520 ka) was 0.45 ‰ lower than that of the late Pleistocene. For the period 1240 to 1320 ka at Site U1385, Birner et al. (2016) suggested that the threshold was 3.2 ‰ on the Cibicidoides scale, which is equivalent to 3.84 ‰ on the Uvigerina scale. This is the same value we have estimated for the entire 1.5-million-year interval, suggesting that the benthic δ¹⁸O threshold has not changed significantly at Site U1385. The existence of an upper δ¹⁸O threshold above which millennial variability is suppressed during peak glacial conditions is less clear from the probability density estimates (Fig. 17). However, several of the late Pleistocene glacial intervals (MIS 2, 4, 6, 10, 12, 16) show a pattern of strong MCV associated with glacial inception that decreases towards full glacial conditions (>4.8 ‰ on the Uvigerina scale), suggesting reduced MCV during maximal glacial conditions.

The physical significance of the benthic δ¹⁸O thresholds that define the millennial window is uncertain. Although several studies have suggested that MCV is related to ice volume, it is not certain which part of the climate–cryosphere system was responsible. Several processes have been suggested to trigger increased MCV, including sea level dropping below a critical sill depth (e.g. Bering Sea; DeBoer and Nof, 2004), the effect of ice sheet height on winds (Zhang et al., 2014), iceberg calving and freshwater fluxes to the oceans, or direct orbital forcing (Friedrich et al., 2010; Zhang et al., 2021; Yin et al., 2021).

McManus et al. (1999) suggested that MCV was enhanced with a sea level lowering of as little as 30 m below modern. This may correspond to a critical sill depth such as the Bering Sea, which has a sill depth of ∼45 m. DeBoer and Nof (2004) proposed that the onset and cessation of flow through the Bering Strait was responsible for the switch between stable and unstable states of glacial versus interglacial climate. A restricted Bering Strait increases the sensitivity of AMOC to freshwater perturbation by blocking the escape route of freshwater to the Pacific via the Arctic (Poppelmeier et al., 2021; Hu et al., 2012a, b). Freshwater can accumulate in the North Atlantic more readily with a closed Bering Strait, thereby increasing surface stratification and leading to AMOC instability.

Because benthic δ¹⁸O also depends on bottom temperature, the threshold could also be related to surface temperature conditions in the source area of deep-water formation. For example, the benthic δ¹⁸O threshold could correspond to crossing the freezing point of seawater in deep-water source areas in the North Atlantic, which would result in increased sea ice formation and shift the region of deep-water formation to the south, where the AMOC is more susceptible to oscillation (Sevellec and Fedorov, 2015).

Galbraith and de Lavergne (2018) suggested that D–O-like variability in AMOC was more likely to occur under a “sweet spot” of interrelated conditions that included low obliquity, low CO₂, and a low-elevation Laurentide ice sheet. By analysing dust flux from the Dome Fuji ice core (Antarctica) over the last 720 ky, Kawamura et al. (2017) also concluded that MCV was more likely during times of intermediate glacial states. Because glacial climate state is ultimately affected by orbital geometry, an inherent link must exist between climate variability on orbital and suborbital timescales (see discussion in Sect. 4.4).

4.6 MCV across the Middle Pleistocene Transition

The MPT is generally considered to have occurred between ∼1200 and 650 ka, although the exact timing is dependent upon the proxy signal and method used to define the shift in frequency (Berends et al., 2021). The benthic δ¹⁸O record involved an increase in amplitude and a shift in the dominant period of glacial–interglacial cycles from 41 kyr before ∼1200 ka to quasi-100 kyr after 650 ka (Clark et al., 2006). Some studies have suggested that MCV was less frequent during the early Pleistocene and increased across the MPT as the size of Northern Hemisphere ice sheets expanded (Larrasoana et al., 2003; Weirauch et al. , 2008; Bolton et al., 2010). Others have found evidence for equally strong millennial variability in the early Pleistocene as in the late Pleistocene (Raymo et al., 1998; McIntyre et al., 2001; Tzedakis et al., 2015; Grützner and Higgins, 2010; Hodell et al., 2008; Birner et al., 2016). Still, others have suggested that MCV (as represented by IRD events) was more frequent but less intense prior to 650 ka because the climate system spent more time in the millennial window during the early Pleistocene (Hodell et al., 2008; Hodell and Channell, 2016) and rarely exceeded the upper benthic δ¹⁸O threshold before 650 ka.

The planktic δ¹⁸O and Zr–Sr records of Site U1385 clearly demonstrate that MCV was strong during glacial stages both before and after the MPT. The main difference across the MPT is that, whereas MCV persists throughout the glacial periods prior to 1200 ka, it is most prevalent on the transitions both into and out of glacial states (i.e. inceptions and terminations) and during times of sustained intermediate ice volume, such as MIS 3. Beginning at 650 ka (MIS 16) following the MPT, MCV is suppressed during some of the strongest glacial periods associated with the growth of oversized continental ice sheets, which Raymo (1997) refers to as “excess ice”.

A change occurred in the orbital modulation of MCV across the MPT, as expressed in changes in the evolutive spectra of the Taner filter–Hilbert transform of the Zr–Sr and planktic δ¹⁸O (Figs. 15 and 16). Prior to ∼900 ka, the amplitude modulation of MCV was dominated by 41 kyr obliquity. Obliquity continues to modulate the amplitude of MCV from 900 to 600 ka but with an increase in precession and the addition of a possible combination tone (28 kyr) of the 41 kyr cycle.

By ∼600 ka, the power of obliquity declines, and the spectra become more complex with greater modulation at lower frequencies (e.g, 100±20 kyr). Hodell et al. (2008) noted a similar change in the amplitude modulation of the Si–Sr IRD proxy at Site U1308 in the central North Atlantic IRD belt when, at ∼650 ka, the power of the 41 kyr obliquity cycle decreased and quasi-100 kyr power increased. Hodell and Channell (2016) also observed that millennial variability in the Si–Sr IRD proxy was proportional to the power in the precessional band, suggesting an amplitude modulation of MCV by precession. Precession plays a greater role in modulating the amplitude of MCV in the late Pleistocene, in agreement with its steady increase in importance throughout the Quaternary (Liautaud et al., 2020).

At 0.65 Ma, the development of massive ice sheets on North America (Batchelor et al., 2019) introduced a new type of MCV related to dynamic instability of the Laurentide ice sheet in the region of Hudson Strait, which was expressed by the occurrence of Heinrich layers in North Atlantic sediment beginning in MIS 16 (Hodell et al., 2008; Hodell and Channell, 2016). Heinrich events tend to occur late in the glacial cycle and are associated with glaciations of long duration (Hodell et al., 2008). They are distinct from background IRD events in their magnitude, frequency, and duration, and their impact on the global climate system was more widespread (Marshall and Koutnik, 2006). MCVs associated with late Pleistocene terminations after 0.65 Ma are closely related to freshwater fluxes from the decay of oversized ice sheets, which play an important role in the progression of glacial terminations (Wolff et al., 2009; Barker and Knorr, 2021).

4.7 Influence of MCV on glacial–interglacial cycles

Ice dynamics may be an effective mechanism for propagating high-frequency variability to longer, orbital timescales (Verbitsky et al., 2019). For example, Siddall et al. (2008) suggested that sustained ice sheet growth through a glacial cycle requires the absence of MCV. Niu et al. (2019) proposed that the presence of strong MCV may prevent an ice sheet from reaching its maximum size owing to surface mass balance effects. If true, then sustained MCV throughout the glacial periods of the early Pleistocene may have prevented ice sheets from growing as large as their late Pleistocene counterparts. In contrast, strong MCV on glacial inceptions may have significantly slowed ice sheet development, giving rise to the sawtooth shape of the late-Pleistocene benthic δ¹⁸O signal. Ice sheets could only reach their maximum size during the latter part of the glacial cycle once MCV was suppressed.

The exact cause–effect relationship between MCV and ice sheet size is difficult to ascertain. Did ice sheets grow larger in the late Pleistocene because MCV was suppressed, or did large ice sheets lead to a suppression of MCV during full glacial conditions? In either case, orbital and millennial-scale variability cannot be considered separately from one another because they interact. Verbitsky et al. (2019) demonstrated that ice sheet non-linearity allows MCV to propagate upscale and influence ice age dynamics. In addition, non-linear ice flow dynamics can propagate downscale and affect the millennial part of the spectrum.

MCV constitutes a source of high-frequency variability on orbital timescales, which may enhance the phase locking of the response of the climate system to orbital forcing (Hodell and Channell, 2016). The theory of stochastic resonance has long been considered as a possible mechanism to explain how the climate system can be synchronized with relatively weak orbital forcing (Benzi et al., 1982). The “noise” for stochastic resonance is often assumed to be random and white. Although MCV is not noise, it provides a source of high-frequency variability in the climate system whose amplitude varies with climate background state – i.e. it is relatively “active” during glacials and “quiet” during interglacials. Such oscillations in amplitude may be relevant for stochastic or coherence resonance in which the signal–noise resonance is important for phase locking (Liu et al., 2018). For example, glacial–interglacial variations during the early Pleistocene may consist of a resonant system in which the intensity of millennial variability responds to obliquity-controlled changes in climate background state, and in turn, changes in the amplitude of MCV may aid in phase locking the climate system to the obliquity period. Stochastic forcing by millennial and centennial climate variability may have also been a crucial factor for the frequency band change associated with the MPT (Mukhin et al., 2019).

4.8 MCV and atmospheric CO₂ variations

Because nearly all the rapidly exchangeable carbon in the ocean–atmosphere system is contained in the deep ocean, atmospheric greenhouse gas variations in ice cores are intimately linked to carbon storage in the deep ocean. Variations in benthic carbon isotopes at Site U1385 demonstrate that the millennial changes in planktic δ¹⁸O are not only a feature of surface climate on the Iberian margin but are crucially linked with changes in deep-water ventilation. Decreases in benthic δ¹³C are associated with increases in planktic δ¹⁸O, indicating reduced ventilation of the deep North Atlantic during cold stadial events. A relationship between atmospheric CO₂ and centennial–millennial events in the North Atlantic exists for the last glaciation and deglaciation (Marcott et al., 2014; Bauska et al., 2021), as well as for older periods such as MIS 6 (Shin et al., 2020) and the MIS 11–10 transition (Nehrbass-Ahles et al., 2020).

We suggest that MCV may be involved in setting the minimum CO₂ values attained during glacial periods. Millennial variability in AMOC provides a mechanism by which deep-sea CO₂ can be degassed to the atmosphere. When MCV was strong during MIS 3, CO₂ varied between 200 and 220 ppm, and the lowest sustained CO₂ levels of 180–190 ppm were only achieved during MIS 2, when MCV was suppressed during peak glacial conditions. By analogy with MIS 3, the persistently strong MCV that occurred throughout the glaciations of the early Pleistocene (1.45–1.25 Ma) may have prevented CO₂ from reaching values as low as those attained during the late Pleistocene because CO₂ was episodically released from the deep-sea reservoir during strong millennial-scale AMOC events. In the early Pleistocene, boron isotope reconstructions suggest that fluctuations in CO₂ varied in phase with obliquity and benthic δ¹⁸O (Chalk et al., 2017; Dyez et al., 2018). The threshold-type behaviour of MCV during the 41 kyr cycles of the early Pleistocene may have served as an important mechanism for linking internal climate dynamics with external astronomical forcing by regulating carbon storage in the deep sea.

Evidence from Site U1385 for an active oceanic thermal bipolar seesaw during most of the prominent stadials during glacials of the 41 kyr world (Birner et al., 2016) supports a similar mechanism of CO₂ degassing via the Southern Ocean as that in MIS 3. Although CO₂ records are fragmentary before 800 ka, there is evidence for elevated glacial CO₂, with minimum values of 220 ppm during glacial periods between 1 and 1.25 Ma during the early MPT (Yan et al., 2019; Higgins et al., 2015; Chalk et al., 2017; Hönisch et al., 2012), and glacial CO₂ values may have been higher still before 1.25 Ma (Yan et al., 2019; Martínez-Botí et al., 2015).

We have emphasized the role that MCV may play in setting atmospheric CO₂ concentrations, but others have suggested that, in contrast, atmospheric CO₂ may have a controlling influence on millennial-scale climate oscillations (Zhang et al., 2017; Vettoretti et al., 2022). Using an Earth system model, Vettoretti et al. (2022) demonstrated that non-linear self-sustained climate oscillations appear spontaneously within an intermediate window of glacial-level atmospheric CO₂ concentrations between ∼190 and 225 ppm.

5 Conclusion

The recognition of MCV in Greenland ice cores in the early 1980s ushered in the study of paleoceanographic records at a resolution that is at least 10 times greater than previous orbital-scale studies. Although the initial focus was on the last deglaciation and MIS 3, several long records of MCV are beginning to emerge (Hodell et al., 2008, 2015; Hodell and Channell, 2016; Barker et al., 2021, 2022), thereby providing an opportunity to document the long-term relationships of climate variability on orbital and millennial timescales and their interactions. Consistent with previous findings, the U1385 record demonstrates that MCV was a persistent feature of intermediate glacial climate states for the last 1.45 Ma, including the 41 kyr world of the early Pleistocene prior to the MPT.

During glacial periods from 1.45 to 1.25 Ma, the amplitude of MCV was strongly modulated by changes in Earth's obliquity and exhibited threshold behaviour typical of a non-linear system. Beginning at 1.2 Ma at the start of the MPT, MCV becomes more focused on glacial inceptions and terminations and periods of intermediate ice volume. One of the recurrent patterns is that strong MCV almost always occurs at glacial inception and continues through the period of ice growth under conditions of declining insolation forced mainly by obliquity and secondarily by precession and CO₂. During the MPT (1.2–0.65 Ma), obliquity continues to influence MCV but in a non-linear fashion, evidenced by the appearance of combination tones (28 kyr) of the 41 kyr cycle (Figs. 15 and 16) in the power spectrum of MCV amplitude modulation. Near the end of the MPT at 650 ka, MCV amplitude modulation by obliquity wanes as quasi-periodic 100 kyr and precession power increases. Precession plays a greater role in modulating the amplitude of MCV in the late Pleistocene, consistent with the steady increase in precession power throughout the Quaternary (Liautaud et al., 2020).

Dansgaard–Oeschger events during MIS 3 are the archetypal example of millennial variability, and considerable effort has been directed towards documenting these events globally, including the use of numerical models to understand their cause(s). MIS 3 is exceptional relative to the other latest Pleistocene glaciations in terms of the high number of millennial events, and there appears to be no period like it during the past 1200 kyr. The strong, continuous millennial variability exhibited during MIS 3 is more similar to the millennial variability observed during glacial cycles of the early Pleistocene from 1440 to 1200 ka (Birner et al., 2016). This similarity is not entirely unexpected considering that benthic δ¹⁸O values were about the same during early-Pleistocene glacial stages as those during MIS 3, indicating that the climate system spent a prolonged time in an intermediate glacial state. Our analysis of MCV at Site U1385 supports the concept of a millennial window or sweet spot defined by a lower benthic δ¹⁸O threshold of ∼2.9 ‰, below which MCV is suppressed during full interglacial conditions. The upper benthic δ¹⁸O threshold is less robust despite the fact that some glacial cycles in the late Pleistocene show a clear pattern of reduced amplitude of MCV as the glacial maximum is approached. Although the exact physical significance of the benthic δ¹⁸O threshold remains uncertain with many candidates (ice volume, ice height, sea level, sea ice, deep-water temperatures, etc.), MCV is strongest during intermediate glacial states.

Climate variability on orbital and suborbital timescales is coupled and interacts in both directions. An example of downscale interaction is the modulation of the amplitude and/or frequency of MCV by Earth's orbital configuration, either through the direct effects of insolation or more indirectly through ice sheet growth. Some MCV may also be related to harmonics or combination tones of the orbital cycle (Hagelberg et al., 1994). MCV can exert an upscale influence on orbital timescales through its effect on ice sheet dynamics (Verbitsky et al., 2019) or on atmospheric CO₂ by changing carbon storage in the deep sea. MCV is also a source of noise on glacial–interglacial timescales that may affect the resonance of internal climate change with external orbital forcing.

In addition to documenting MCV, the planktic and benthic isotope records from Site U1385 provide unprecedented detail of the amplitude and shapes (waveforms) of the glacial cycles on orbital timescales for the last 1.45 Ma. We emphasize that our record is from a single site and should be compared with other records from higher latitudes in the North Atlantic (e.g. Barker et al., 2021, 2022) and elsewhere (Sun et al., 2021) in order to map geographical differences over time and to develop confidence in the palaeoceanographic interpretations set out here. This study is also limited to the last 1.45 Ma, and we cannot determine the extent to which MCV was present during glacial periods beyond this time. In late 2022 (12 October–12 December), International Ocean Discovery Program Expedition 397 returned to the Iberian margin and extended the record of Site U1385 to 4.5 Ma in the early Pliocene (Hodell et al., 2023). The new sediment cores recovered during Expedition 397 (Iberian margin paleoclimate) will document how orbital and millennial variability co-evolved as climate background state changed from the warm conditions of the early and middle Pliocene through the intensification of Northern Hemisphere glaciation during the late Pliocene and Quaternary. Understanding the interactions of climate on orbital and suborbital timescales will lead to a fuller understanding of the mechanisms responsible for the Quaternary ice ages.