Millennial-scale atmospheric CO₂ variations during the Marine Isotope Stage 6 period (190–135 ka)

Using new and previously published CO₂ data from the EPICA Dome C ice core (EDC), we reconstruct a new high-resolution record of atmospheric CO₂ during Marine Isotope Stage (MIS) 6 (190 to 135 ka)  the penultimate glacial period. Similar to the last glacial cycle, where high-resolution data already exists, our record shows that during longer North Atlantic (NA) stadials, millennial CO₂ variations during MIS 6 are clearly coincident with the bipolar seesaw signal in the Antarctic temperature record. However, during one short stadial in the NA, atmospheric CO₂ variation is small (∼5 ppm) and the relationship between temperature variations in EDC and atmospheric CO₂ is unclear. The magnitude of CO₂ increase during Carbon Dioxide Maxima (CDM) is closely related to the NA stadial duration in both MIS 6 and MIS 3 (60–27 ka). This observation implies that during the last two glacials the overall bipolar seesaw coupling of climate and atmospheric CO₂ operated similarly. In addition, similar to the last glacial period, CDM during the earliest MIS 6 show different lags with respect to the corresponding abrupt CH₄ rises, the latter reflecting rapid warming in the Northern Hemisphere (NH). During MIS 6i at around 181.5±0.3 ka, CDM 6i lags the abrupt warming in the NH by only 240±320 years. However, during CDM 6iv (171.1±0.2 ka) and CDM 6iii (175.4±0.4 ka) the lag is much longer: 1290±540 years on average. We speculate that the size of this lag may be related to a larger expansion of carbon-rich, southern-sourced waters into the Northern Hemisphere in MIS 6, providing a larger carbon reservoir that requires more time to be depleted.

1 Introduction

Ice core studies allow us to considerably extend our knowledge about natural climate–carbon cycle feedbacks by directly reconstructing atmospheric CO₂ from gas preserved in Antarctic ice sheets (Lüthi et al., 2008; Petit et al., 1999). Comparing atmospheric CO₂ records from Antarctic ice cores with proxies of paleoclimate helps us to understand how atmospheric CO₂ was controlled by carbon exchange with the ocean and land reservoirs on orbital to centennial timescales (Ahn and Brook, 2008, 2014; Bereiter et al., 2012; Higgins et al., 2015; Lüthi et al., 2008; Marcott et al., 2014; Petit et al., 1999).

Previous work on polar ice core records revealed that temperature variations in Greenland and Antarctica on millennial timescales appear to be a pervasive feature during the last glacial period. While Antarctic temperature varied gradually, Greenland temperature changes occurred abruptly. A phase difference can be observed between millennial-scale variations of temperature in the NH and SH (Northern Hemisphere and Southern Hemisphere, respectively), which is referred to as the bipolar seesaw phenomenon (Blunier and Brook, 2001; Pedro et al., 2018; Stocker and Johnsen, 2003). Potential triggers for this climatic variability on the millennial scale are fresh water perturbation in the North Atlantic (NA) or alterations of sea ice extent, surface temperature and salinity in the NA (Bond et al., 1992; Broecker et al., 1992; Heinrich, 1988; McManus et al., 1998), which may reduce the strength of the Atlantic Meridional Overturning Circulation (AMOC). This would cause a reduction in heat transport from the SH to the NH, which leads to an abrupt cooling in the NA region and a gradual warming in the SH (Stocker and Johnsen, 2003) and the opposite behavior when AMOC is strengthened.

Existing CO₂ records show the presence of millennial-scale oscillations on the order of 20 ppm over the last glacial period (Ahn and Brook, 2008, 2014; Bereiter et al., 2012), which generally co-vary with the major water isotope (δD) variations in Antarctic ice cores reflecting Antarctic temperature variations (Ahn and Brook, 2008; Bereiter et al., 2012) (Fig. 1). During cold periods in the NA, referred to as NA stadials, atmospheric CO₂ increased continuously and in parallel to Antarctic temperature increase. Again, at the onset of warming in Greenland, atmospheric CO₂ started to decrease (Ahn and Brook, 2008; Bereiter et al., 2012), generally in line with a co-occurring, slow Antarctic cooling. However, the CO₂ decrease did not always start at exactly the same time as the onset of the abrupt warming in the NA, and the lag itself varied. For example, during Marine Isotope Stage (MIS) 3, Carbon Dioxide Maxima (CDM) lagged behind abrupt temperature change in Greenland by 870±90 years. During MIS 5, the lag of CDM with respect to abrupt temperature warming in the NH was only about 250±190 years (Bereiter et al., 2012). Atmospheric CO₂ variations on millennial scales are thought to be related to the role of the Southern Ocean in carbon uptake and deep-ocean ventilation on millennial timescales (Fischer et al., 2010; Marcott et al., 2014; Schmitt et al., 2012; Sigman and Boyle, 2000; Toggweiler et al., 2006). In addition, atmospheric CO₂ can be affected by changes in the AMOC, which affects the ventilation of carbon from the deep ocean (Denton et al., 2010; Sigman et al., 2007). However, the mechanisms responsible for these oscillations are still a matter of debate (Bouttes et al., 2012; Bozbiyik et al., 2011; Gottschalk et al., 2019; Menviel et al., 2014).

Comparing CO₂ changes on millennial timescales during the past two glacial periods, MIS 3 (60–27 ka) and early MIS 6 (185–160 ka), can provide us with a better understanding of the carbon cycle, due to the similarities but also differences of climate conditions and events during the last two glacial periods (Fig. 1). Proxy evidence indicates that the states of several important components of the climate–carbon cycle were not the same between MIS 3 and MIS 6. Sea ice cover in the South Atlantic was more extensive in MIS 6, and sea surface temperature in the South Atlantic is thought to have been lower (Gottschalk et al., 2020). The bipolar seesaw phenomenon has also been observed during the early MIS 6 period (Cheng et al., 2016; Jouzel et al., 2007; Margari et al., 2010). However, the bipolar seesaw events during MIS 6 are longer than those found during MIS 3. Events of massive iceberg discharge into the NA, which are thought to have driven millennial-scale changes in the meridional overturning circulation during MIS 3 (de Abreu et al., 2003; McManus et al., 1999), appear to be much more frequent during MIS 3 than during MIS 6. During the early MIS 6, iceberg discharge was muted, and during the time period around 175 ka summer insolation levels in the NH approached interglacial values (Berger, 1978). Due to the stronger NH summer insolation, the Intertropical Convergence Zone (ITCZ) had shifted to the north, which intensified monsoon systems in low-latitude regions, such as in Asia, the Apennine Peninsula and the Levant (Ayalon et al., 2002; Bard et al., 2002; Cheng et al., 2016). This may have led to a weaker overturning circulation due to the reduction of the density of the North Atlantic surface water, making the AMOC cell shallower during MIS 6 than during MIS 3 (Gottschalk et al., 2020; Margari et al., 2010).

Figure 1 Proxy data over the last 250 ka. (a) Ice-rafted debris (IRD) input in the Iberian Margin core MD95-2040 (de Abreu et al., 2003). (b) The 21 June insolation for 65^∘ N (Berger, 1978). (c) The δ¹⁸O_calcite from Sanbao cave, indicative of the strength of the East Asian monsoon (Cheng et al., 2016). (d) Dust flux in the EDC ice core (Lambert et al., 2012). (e) Published atmospheric CH₄ in EDC (dark green dots) (Loulergue et al., 2008) and a composite atmospheric CH₄ record from EDC derived in this study (light yellow dots). (f) Composite CO₂ record from the EDC ice core derived in this study (orange dots) and a published composite CO₂ record from Antarctic ice cores (dark blue dots) (Bereiter et al., 2015). (g) The δD in the EDC ice core, Antarctica (Jouzel et al., 2007). Vertical grey bars indicate the timing of Heinrich events during the last glacial period. The event numbers are indicated at the bottom. The marine isotope substage numbers are written at the bottom (Railsback et al., 2015).

In order to investigate whether the different climate conditions between MIS 6 and MIS 3 could have impacted the relationship between atmospheric CO₂ and climate, we measured 177 new data points of atmospheric CO₂ concentrations from the EPICA Dome C (EDC) ice core (75^∘06^′ S, 123^∘24^′ E) spanning MIS 6 (190 to 135 ka). This new CO₂ data measured within the scope of this study is complemented with two published datasets (Lourantou et al., 2010; Schneider et al., 2013) and one unpublished CO₂ dataset measured in 2003. We use all of the available data from the EDC ice core to compile a composite dataset of atmospheric CO₂ covering MIS 6. In this work, we significantly improve existing records previously obtained from the Vostok ice core (Petit et al., 1999). The new composite dataset from the EDC core provides a temporal resolution of ∼230 years during MIS 6, as compared to 1000 years in the Vostok record. We also improved the relative age uncertainties between ice and gas in the EDC core using δ¹⁵N-based estimates of firn column thickness to better constrain leads and lags between peaks in δD (Antarctic Isotope Maxima, or AIM) in the EDC ice core and atmospheric CO₂ concentrations during the early MIS 6, thus establishing a new chronology using new δ¹⁵N data during the early MIS 6 in this study and published data from Landais et al. (2013). Finally, we improve the temporal resolution of existing CH₄ data from EDC (Loulergue et al., 2008) from 600 to ∼350 years to be able to more precisely calculate the shift of CDM relative to the rapid climate change in the NH. To avoid the age uncertainties between proxy data and atmospheric CO₂ data, CH₄ measurements are used as a time marker of rapid warming in the NH, as over the last glacial period CH₄ and Greenland temperatures are assumed to be synchronous, with a lag of CH₄ of less than ∼50 years on average (Baumgartner et al., 2014; Rosen et al., 2014).

2 Methods

2.1 CO₂ measurements

In this study, we measured atmospheric CO₂ from 177 depth intervals of the EDC ice core at the Institut des Géosciences de l'Environnement (IGE), France and Climate and Environmental Physics (CEP), Physics Institute, University of Bern, Switzerland. The majority of all atmospheric CO₂ samples (i.e. 150) were measured using the ball mill dry-extraction system coupled to a gas chromatograph at IGE (Schaefer et al., 2011). Each ball mill data point presented in this study corresponds to a single 40 g ice sample, which was measured five times by gas chromatography (five consecutive injections of the same extracted gas). Approximately 5 mm of ice was trimmed from the ice core surfaces before extraction in order to remove the external part that could be potentially contaminated with drilling fluid or might have been subject to gas loss during storage in the freezer (Bereiter et al., 2009). The CO₂ measurements were referenced to a secondary gas standard (synthetic air from Air Liquide, Alphagaz 28416000) containing 233.7±0.4 ppm of CO₂ in dry air, which was referenced to two primary standards (238.34±0.04 from NOAA, CB09707, and 260.26±0.2 from CSIRO, CSIRO1677).

Blank tests using 40 g of artificial bubble-free ice were conducted every 10 measurements to quantify the precision of the system and to correct for the CO₂ contamination caused by the crushing process. Blank tests were conducted in two steps: first, to validate the baseline of the system, a gas standard with 233.7±0.4 ppm was injected over the bubble-free ice in the cell. The gas was then left to equilibrate in the cell for 10 min. The gas was analyzed twice by successive injections into the extraction line and sample loop. Afterwards, the bubble-free ice was crushed and the gas was analyzed five more times. The difference between the results before and after crushing was considered the contamination effect caused by the crushing process. These values were used to estimate the precision of the system. Measured CO₂ should be corrected for contamination caused by the analytical procedure by comparing measured CO₂ of the blank tests with the standard gas value. However, it is not feasible to correct CO₂ concentrations directly. The CO₂ mole fraction is calculated as the ratio between partial pressure of CO₂ and total pressure in the measurement line, which is a relative value. Thus, CO₂ concentration and the concentration of CO₂ contamination are dependent on the total pressure. To avoid this dependence, we corrected the absolute value (partial pressure) of CO₂ in the air by the expected partial pressure of CO₂ contamination, as estimated from the blank tests.

For this study, we used four different extraction chambers to hold the ice core samples during crushing and measurement. Each chamber showed different contamination levels. Therefore, blank tests were conducted on each chamber. The data were corrected by the average of each chamber. For these measurements, the precision of the system was estimated to be ∼1 ppm on average. On average, the extraction effect correction corresponds to a reduction of the measured CO₂ concentration by 1.7±1.0 ppm (1σ).

Additionally, we measured 65 samples from 27 depth intervals using the Centrifugal Ice Microtome (CIM) (Bereiter et al., 2013) at CEP. At each depth interval replicates were made and analyzed (2.4 on average). Approximately 5 mm of the ice surface were trimmed, and ∼8 g of samples was analyzed on average. The pooled standard deviation of replicates was ∼1 ppm. The WMOX2007 mole fraction scale (Tans et al., 2017) was used as a reference for the CO₂ measurements, via four different primary dry air standards for atmospheric CO₂, which were calibrated at the NOAA Earth System Research Laboratory. The standard sample concentrations varied between 192.44 and 363.08 ppm (Nehrbass-Ahles et al., 2020).

Each CO₂ record was corrected for gravitational fractionation, using the δ¹⁵N isotope ratio (Craig et al., 1988). To this end, 88 new data points together with existing δ¹⁵ (Landais et al., 2013) covering the late MIS 6 (156.4–139.2 ka) were used. δ¹⁵N data were linearly interpolated in age to each corresponding CO₂ data point. On average, the correction corresponds to removing 1.2±0.1 ppm from the measured CO₂. In total, 177 individual ice samples were measured for CO₂ in the depth range from 2036.7 to 1787.5 m along the EDC ice core, corresponding to the time period from 189.4 to 135.4 ka (Bazin et al., 2013).

2.2 CH₄ measurements

We measured the atmospheric CH₄ content of 63 ice core samples, using the wet extraction method at IGE as described in detail in Spahni et al. (2005). This allowed us to improve the temporal resolution of existing CH₄ data (Loulergue et al., 2008) from ∼600 to 360 years on the AICC2012 chronology (see Fig. S1 in the Supplement). The precision of the system is ∼11 ppb on average, which was estimated by blank tests.

The previous CH₄ dataset (Loulergue et al., 2008) from EDC has been produced at both IGE and CEP. CH₄ measured at CEP used to be systematically higher than CH₄ measured at IGE by 6 ppb (Loulergue et al., 2008). The offsets are due to differences in corrections for contamination caused by the analytical procedure between the datasets. Previously, 6 ppb had been added to the data obtained at IGE (Loulergue et al., 2008).

In this study, a new data correction was applied to data measured at IGE, and it is no longer necessary to add 6 ppb to the new data. The systematic offset between the new data and the previously published data points (Loulergue et al., 2008) is only 1.7±2.4 ppb (n=63, the standard error of the mean) during MIS 6, (the datasets agree with each other within the measurement uncertainty).

2.3 Nitrogen isotopes

Isotopes of molecular nitrogen in air bubbles were measured by a melting technique at the Laboratoire des Sciences du Climat et de l'Environnement (LSCE), France. The gas was extracted from the ice by a wet extraction technique, and the released air was analyzed by a dual-inlet mass spectrometer (Delta V Plus; Thermo Scientific). The analytical method and data correction are described in detail in Bréant et al. (2019). In total, 151 samples from 88 depth intervals (63 duplicates) between the depths of 2124.7 and 1875.0 m below the surface were measured, corresponding to the time interval from 205 to 154 ka (Fig. S2 in the Supplement). The average resolution on the AICC2012 chronology is ∼580 years.

2.4 Gas age revision by estimating Δdepth from δ¹⁵N

The water isotopic signature (δD), is, unlike CO₂, measured on the ice matrix. Air enclosed in an ice core moves through the porous firn layers at the top of the ice sheet by molecular diffusion and becomes gradually trapped in the ice at the so-called Lock-In Depth (LID, around 100 m below the surface in the case of the EDC ice core). An age difference, thus, exists between the ice and the air at a given depth. For the conditions of the EDC ice core, the age difference can reach up to 5 kyr during glacial maxima and is associated with a large uncertainty (several hundred years; Loulergue et al., 2007). The δ¹⁵N of molecular nitrogen in air bubbles can be used to determine the LID and to calculate the depth difference between synchronous events in the ice matrix and air bubbles, called delta depth (Δdepth), thus creating a more precise relative chronology of the gas – with respect to the ice-phase (Parrenin et al., 2013). We use the Δdepth calculation to adjust the gas chronology of EDC, while the ice chronology from AICC2012 remained unchanged.

We calculate the height of the firn column, h, from δ¹⁵N of N₂ measurements (Craig et al., 1988; Dreyfus et al., 2010; Sowers et al., 1989) using the following equation:

$$\begin{array}{}\text{(1)}& h=h_{\mathrm{conv}}+({\mathit{\delta }}^{\mathrm{15}}\mathrm{N}-\mathrm{\Omega }(T\left)\mathrm{\Delta }{T}_{\mathrm{diff}}\right){\left({\displaystyle \frac{\mathrm{\Delta }m\cdot g\cdot \mathrm{1000}}{RT}}\right)}^{-\mathrm{1}}.\end{array}$$

In this equation, ΔT_diff is the temperature difference between the top and the bottom of the diffusive zone as estimated by the Goujon–Arnaud model, where surface temperature and accumulation are estimated from the stable water isotope record (Loulergue et al., 2007). Ω (T) is the thermal diffusion sensitivity, which has been estimated from laboratory measurements by Grachev and Severinghaus (2003). Δm is the mass difference between ¹⁴N and ¹⁵N (kg mol⁻¹), g is the gravitational acceleration (9825 m s⁻² for Antarctica) (Parrenin et al., 2013), and R is the universal gas constant (8314 J mol⁻¹K⁻¹). Finally, h_conv is the height of the convective zone at the top of the firn column, which is considered negligible at EDC according to current observations (Landais et al., 2006). The variation of the convective height is related to changes in wind stress. According to Krinner et al. (2000), wind on the East Antarctic plateau varied little during the LGM (Last Glacial Maximum), and we assume that this is also the case for late MIS 6. Moreover, the convective zone was confirmed to be very thin during the last deglaciation by Parrenin et al. (2012). Thus, we assume that h_conv is negligible during MIS 6. Δdepth is calculated from the height of the air column using a constant average firn density and a modeled vertical thinning function suggested by Parrenin et al. (2013).

Raw δ¹⁵N data cannot be used directly to calculate Δdepth because bubbles at a given depth in the ice sheet can be trapped at slightly different times, due to small physical variations that affect the gas trapping process. This can lead to age scale inversions on depth intervals up to 1 cm (Etheridge et al., 1992; Fourteau et al., 2017). These layer inversions should not strongly affect our record, as its resolution is larger than the scale of age inversion events (Fourteau et al., 2017). However, the change in Δdepth between two different depths (z) in the ice core, denoted $\partial \mathrm{\Delta }\text{depth}/\partial z$, deduced from the raw δ¹⁵N data, shows some values higher than 1, that could correspond to small stratigraphic inversions in the gas phase (Fig. 2). Therefore, a three-point moving average weighted by the time difference between a point and its two neighbors was applied to the δ¹⁵N dataset. The weights for the three points are equal if the time difference is less than or equal to 500 years, which is close to the average sampling resolution of the δ¹⁵N dataset. Otherwise, the neighboring data points are weighted by 500∕ΔT, where ΔT is the time difference. This avoids assigning too much weight to neighboring points where the resolution in the record is lower, which would lead to local over-smoothing.

Figure 2 (a) $\partial \mathrm{\Delta }\text{depth}/\partial z$ as a function of depth. Red dots represent raw δ¹⁵N measurements, and blue dots represent a weighted three-point running mean (see text for details). The dashed vertical grey line indicates $\partial \mathrm{\Delta }\text{depth}/\partial z=\mathrm{1}$. (b) Δdepth (bold line) for the EDC ice core from 1787.5 to 1870.2 m below the surface, deduced from δ¹⁵N and the thinning function calculated in this depth range. The two dashed lines correspond to the analytical uncertainties.

Δdepth as estimated using the three-point moving average weighted by the time difference is shown in Fig. 2. The difference of the Lock-In Depth in Ice Equivalent (LIDIE) calculated on the AICC2012 age scale (Bazin et al., 2013) and deduced from δ¹⁵N in this study is 0.5±3.0 m, or 0.7±5.6 % on average (see Fig. S3 in the Supplement). The AICC2012 LIDIE was calculated using a background scenario derived from δD, using a linear relationship between δD and δ¹⁵N. Our results show this relationship to be relatively unbiased but not entirely exact (Fig. S3 in the Supplement). The Δdepth values are used to update the EDC gas chronology from the original AICC2012 age scale, while the ice chronology remains unchanged. For MIS 6, this method significantly reduces the relative age uncertainty between air and ice to 900 years on average with respect to the original AICC2012 chronology (see Fig. S4 in the Supplement).

2.5 Age scale of the MD01-2444 marine sediment record

The MD01-2444 marine sediment core from the Iberian Margin (Margari et al., 2010) represents an important archive for the interpretation of our CO₂ record, as it is well resolved during MIS 6, and it provides high-resolution benthic and planktonic foraminiferal records. The original MD01-2444 age scale was built by matching the benthic δ¹⁸O record (Shackleton et al., 2000) to the δD record from EDC on the EDC03 ice age scale (Parrenin et al., 2007) using 11 tie points (Margari et al., 2010). In order to make optimal use of the record, its relative chronology with respect to EDC requires additional improvement using the latest version of the EDC age scale. Thus, in this study, we recalculated the age of the tie points using the AICC2012 chronology (Bazin et al., 2013), and the age of the sediment record was linearly interpolated between the tie points (Table S1 in the Supplement).

2.6 Definition of NA stadial duration

Due to the absence of a Greenland temperature record for MIS 6, the durations of the six NA stadials must be identified using other proxies. Here, we used two methods using both ice core and sediment core proxies. First, we estimated the durations of the six NA stadials using δ¹⁸O of planktonic foraminifera and tree pollen in MD01-2444, which reflect temperature variability in the NH (Margari et al., 2010) (Fig. 3). The midpoint of the stadial transitions in both δ¹⁸O of planktonic foraminifera and tree pollen in MD01-2444 were used to identify the NH stadial transitions. The time interval between two stadial transition points were defined as the NA stadial duration. Using this approach the duration of stadial at 6ii in the NH could not unambiguously be determined due to small variations in the foraminifera and the pollen record. However, the average age difference between the durations identified using the two records is only 205 years, which is less than the sampling resolution of MD01-2444 during MIS 6. The uncertainty of the timing of each stadial transition was estimated as half of the temporal difference between maxima and minima of δ¹⁸O of planktonic foraminifera. The mean values for stadial onsets and terminations are shown as dashed lines in Fig. 3. However, not all of the stadial durations during MIS 6c to 6d are entirely unambiguous using this method.

Figure 3 The durations of the six NA stadials during MIS 6: (a) δD of the EDC ice core (Jouzel et al., 2007), (b) synthetic Greenland δ¹⁸O_ice record (Barker et al., 2011), (c) tree pollen percentage in the MD01-2444 core (Margari et al., 2010) and (d) δ¹⁸O of planktonic foraminifera in the MD01-2444 core (Margari et al., 2010). Proxy data shown here are given on the AICC2012 age scale. Red lines indicate the midpoints of the stadial transition of both δ¹⁸O of planktonic foraminifera and tree pollen in MD01-2444. Light green bars indicate the uncertainty of the duration of each stadial transition, estimated as half the temporal difference between maxima and minima of δ¹⁸O of planktonic foraminifera before and after the transition. Red dots indicate minima and maxima of δD in the EDC ice core as selected in this study. The event numbers are indicated at the top.

Second, we recalculated the durations of the six NA stadials during MIS 6 period using a method developed by Margari et al. (2010). Figure 3 also shows a synthetic Greenland δ¹⁸O_ice record (Barker et al., 2011) and Antarctic δD variations on the AICC2012 age scale. The interval between the maximum and the preceding minimum of δD in the EDC record can also be used to estimate the duration of the stadial transitions (Gottschalk et al., 2020; Margari et al., 2010). In most cases, the synthetic Greenland δ¹⁸O_ice record and the interval between the maximum and the preceding minimum of δD in the EDC record confirm the definition of NA stadials selected by δ¹⁸O of planktonic foraminifera and tree pollen in MD01-2444. However, the duration of the NA stadial in MIS 6iii is not clearly confirmed by synthetic Greenland δ¹⁸O_ice and δD in the EDC (Fig. 3).

The interval between each maximum and the preceding minimum of δD in the EDC record was calculated to estimate the stadial durations using the second method (Gottschalk et al., 2020; Kawamura et al., 2017; Margari et al., 2010). Maxima and minima in the δD record were selected by finding zero values in the second (Savitzky–Golay filtered) derivative of the data (the same method we used to pick minima and maxima of atmospheric CO₂ and temperature; Figs. S5–S6, Tables S2–S3 in the Supplement). In one case (the minimum before MIS 6i) two potential minima were selected by this method; these values were then averaged.

The red dots and error bars in the EDC δD record in Fig. 3 show the estimated minima and maxima of temperature corresponding to stadial transitions using this method, along with their uncertainties. However, using this tool, durations of MIS 6ii and 6i are apparently overestimated due to ambiguities concerning the definition of the maximum in MIS 6i and minimum in 6ii.

Neither our method nor that of Margari et al. (2010) can be considered absolutely correct. To account for the differences between the two methods, we took the stadial duration to be the mean of the durations estimated by the two methods (Table S4 in the Supplement).

3 Results

3.1 Data compilation

We measured 177 samples of atmospheric CO₂ from EDC during MIS 6 using two extraction systems (the ball mill at IGE and CIM system at CEP). To improve the resolution of our new dataset even further, we made a composite dataset by aligning previous sets of CO₂ measurements made over the MIS 6 period on the EDC ice core to our new data. First, we compared the two existing CO₂ datasets and the new CO₂ dataset from EDC (Fig. 4 and Table S5 in the Supplement). There are two published CO₂ datasets for EDC during MIS6 – the first measured using the ball mill system at IGE (Lourantou et al., 2010) and the second using a sublimation extraction system at CEP (Schneider et al., 2013). We also added unpublished data measured in 2003 using the CEO needle cracker measurement system in 2003 (Monnin et al., 2004; Siegenthaler et al., 2005) (see Supplement for details). All records are on the AICC2012 air age scale (Bazin et al., 2013). All datasets are corrected for the gravitational fractionation effect using the new δ¹⁵N data in our study. In total, the datasets contain 237 CO₂ measurement points. Two samples were excluded because of system operator error. Figure 4 shows CO₂ concentrations measured by the ball mill system, sublimation (Schneider et al., 2013), the CIM and the needle cracker. Concentrations from the CIM, sublimation and the needle cracker are systematically higher than CO₂ concentrations measured using the ball mill system (Table S5 in the Supplement).

Figure 4 Atmospheric CO₂ from EDC and Vostok ice cores, compared to δD of the ice in the EDC core (temperature proxy) during 190–135 ka. Blue dots represent atmospheric CO₂ from EDC as measured by the ball mill system (this study). The error bars of the new CO₂ data using the ball mill extraction system indicate the standard deviation of five consecutive injections of the gas extracted from each sample into the gas chromatograph added to the precision of the measurement estimated by the reproducibility of the control measurement (∼0.8 ppm) using a quadratic sum. Yellow dots represent atmospheric CO₂ from EDC measured using the ball mill system (Lourantou et al., 2010). The error bars indicate the standard deviation of five consecutive injections of the gas extracted from each sample into the gas chromatograph. Red triangles represent atmospheric CO₂ from EDC measured using the needle cracker (this study). The error bars of data points with replicates indicate the standard deviation of the mean of replicates from the same depth interval. Black inverted triangles represent atmospheric CO₂ from EDC measured using the CIM (this study). The error bars of data points with replicates indicate the standard deviation of the mean of replicates from the same depth interval. Green diamonds represent atmospheric CO₂ from EDC measured using sublimation extraction (Schneider et al., 2013). The error bars of data points with replicates indicate the standard deviation of the mean of replicates from the same depth interval (Schneider et al., 2013). Grey dots represent atmospheric CO₂ from the Vostok ice core (Petit et al., 1999). The error bars of the CO₂ data from Vostok (Petit et al., 1999) show the estimated overall accuracy for CO₂ measurements. The grey line represents δD of ice in the EDC core (Jouzel et al., 2007).

The offsets between each additional dataset and our data were calculated and corrected using a Monte Carlo procedure to account for uncertainty and a Savitzky–Golay filter used to remove noise from our dataset (Supplement). The offsets between the multiple datasets may be at least in part linked to differences in extraction efficiency between the measurement methods. The sublimation and CIM systems have high extraction efficiency on clathrates and should therefore present more unbiased baseline CO₂ values. However, since these datasets are of much lower resolution, we have aligned all datasets to the baseline value of our ball mill dataset.

Despite the systematic offsets, the other datasets confirm the millennial-scale variations shown in the data from this study (Fig. 4 and Fig. S7 in the Supplement). Although none of the individual additional datasets is of high enough resolution to show millennial-scale variations in detail, when aligned to our data the new data follow the millennial-scale variations with very few outliers. The new composite atmospheric CO₂ record from EDC is shown in Fig. 5.

Figure 5 Composite atmospheric CO₂ (left axis) from the EDC ice core (this study) compared to the EDC water isotopic record (right axis) (Jouzel et al., 2007). The blue line indicates a Savitzky–Golay filter using a 150-year cut-off period (dotted blue line). Vertical dotted blue lines indicate the six CDM events that we identify during the early MIS 6. The numbers of the CDM events are indicated at the top of the figure.

3.2 Comparison with Vostok CO₂ record for MIS 6

The new composite atmospheric CO₂ record from EDC was compared to the existing CO₂ data from the Vostok ice core, measured using the ball mill system (Fig. S7 in the Supplement), where the CO₂ record from Vostok was aligned to the AICC2012 gas age scale (Bazin et al., 2013). Atmospheric CO₂ data from Vostok are corrected for the gravitational fractionation effect using the existing δ¹⁵N data (Bender, 2002).

Atmospheric CO₂ variability in the EDC ice core shows similar general patterns as those retrieved from the Vostok ice core, although previously unidentified features are resolved due to the improved temporal resolution. CO₂ concentrations from Vostok appear systematically higher than those from EDC by 4.6±3.0 ppm on average and the difference increases to more than 10 ppm at the beginning of Termination 2 at around 135 ka.

Three main mechanisms may explain the offsets between the EDC and Vostok measurements: firstly, contamination due to the extraction system was assumed to be negligible at the time of the Vostok sample measurements at IGE (Petit et al., 1999). Based on recent ball mill measurements conducted within the scope of this study, we presume that the neglected correction concerning the data published by Petit et al. (1999) amounts to approximately 1.7 ppm, reducing the average apparent offset to about 3 ppm. Secondly, part of this remaining offset between the new EDC and previously published Vostok record may be related to age scale uncertainties due to a limited number of stratigraphic tie points between the two cores (Bazin et al., 2013). In particular, this effect may explain part of the larger offset during the glacial termination in air younger than 140 ka (Fig. S7 in the Supplement). Thirdly, the ball mill system has a different extraction efficiency depending on the presence of bubbles and/or clathrates in the ice sample, which may affect the accuracy of the reconstructed absolute mean CO₂ level. When the air is extracted from an ice core sample where bubbles and clathrates co-exist, different proportions of bubbly and clathrate ice may lead to biased CO₂ concentrations (Lüthi et al., 2010; Schaefer et al., 2011). The Vostok measurements were made on recently drilled ice, in which clathrates had less time to transform into secondary bubbles (Vladimir Lipenkov, 2015 personal communication). Accordingly, at least part of the approximately 3 ppm offset may be due to these systematic extraction differences (see Supplement).

In summary, the individual datasets used for our data compilation and the Vostok ice core record all show the same millennial CO₂ variability over MIS 6. Analytical and ice-relaxation-related offsets could be removed to synthesize all records. This allows us to compile a precise and high-resolution record; however, the true absolute mean level of atmospheric CO₂ is only known to be better than 5 ppm, due to this offset correction. However, this has no impact on the conclusions drawn from the millennial CO₂ variability, which is the focus of this study.

3.3 Atmospheric CO₂ variability on millennial timescale

Margari et al. (2010) suggested that MIS 6 can be divided into three sections depending on the degree of climatic variability observed in δD (indicative of Antarctic climate variability) and CH₄ (reflecting NA climate variability) in EDC: early (185.2–157.7 ka), transition (157.7–151 ka) and late MIS 6 (151–135 ka) (Figs. 5 and 6). Climatic oscillations on millennial timescales are pervasive during the early MIS 6 period (185–160 ka) (Barker et al., 2011; Cheng et al., 2016; Jouzel et al., 2007; Margari et al., 2010, 2014), which is similar to MIS 3 (Figs. 1 and 6). However, during the late MIS 6 period, i.e., the Penultimate Glacial Maximum (PGM), millennial variability is subdued and more resembles the climate variability on millennial timescales during MIS 2 (Figs. 1 and 6). During the transitional period from 157 to 151 ka, δD in EDC slowly increased (Jouzel et al., 2007). Like δD in EDC, CO₂ variations on millennial timescales are pervasive during the early MIS 6 period (185–157.7 ka). During the transitional period from 157 to 151 ka, atmospheric CO₂ increased slowly, while during the late MIS 6 period CO₂ variation is subdued (Fig. 5).

Figure 6 Comparison of climate proxies with atmospheric CO₂ during MIS 6 period. Vertical dotted blue lines indicate the six CDM events during the early MIS 6: (a) 21 June insolation at 65^∘ N (Berger, 1978), (b) ice-rafted debris (IRD) input in the Iberian Margin core MD95-2040 (de Abreu et al., 2003), (c) atmospheric CH₄ in the EDC ice core (Loulergue et al., 2008; this study), (d) δ¹⁸O of planktonic foraminifera in the Iberian Margin marine core MD01-2444 (Margari et al., 2010), (e) δ¹⁸O of benthic foraminifera in the Iberian Margin marine core MD01-2444 (Margari et al., 2010), (f) δD of the EDC ice core (Jouzel et al., 2007) and (g) our new composite CO₂ record during the MIS 6 period. The numbers of CDM events are indicated at the top.

While there was an indication of millennial CO₂ variability in the Vostok record, the resolution of that record did not allow us to constrain this variability. We now present clear evidence of millennial variability of CO₂ concentrations during MIS 6 that are associated with Antarctic Isotope Maxima (AIM) events (EPICA Community Members, 2006), thanks to the improved time resolution and precision of the obtained CO₂ data (Fig. 5). To better discuss millennial variability in MIS 6, a Savitzky–Golay filter with a 150-year cut-off period was used to filter out centennial-scale variability and noise (see Fig. 5 and Fig. S5 in the Supplement). Five prominent and one subdued CO₂ variations were detected in atmospheric CO₂ during early MIS 6 (Fig. 5 and Fig. S5 and Table S2 in the Supplement), where we name CDM according to the numbering by Margari et al. (2010) and Gottschalk et al. (2020) (Fig. 5). The five prominent peaks are observed at 160.7±0.3 (CDM 6vi), 164.2±0.3 (CDM 6v), 169.6±0.2 (CDM 6iv), 174.3±0.2 (CDM 6iii) and 181.3±0.1 ka (CDM 6i). Note that the provided uncertainty was calculated with respect to the position of each maximum and does not include the absolute age uncertainty of the ice core (in each case around 3.0 kyr, 1σ) (Bazin et al., 2013). We can also identify one low-amplitude CO₂ peak at around 150 ka, representing another potential candidate for a CDM (Fig. 5). This atmospheric CO₂ variation is of triangular shape and follows the δD pattern. The change of direction is also associated with a CH₄ peak. This variation has analogues in MIS 4 and MIS 10 (Jouzel et al., 2007; Nehrbass-Ahles et al., 2020).

Each CDM is associated with an AIM event. The short AIM 6ii event corresponds to CDM 6ii at around 178 ka. CDM 6ii has an amplitude of only ∼4 ppm, and the CO₂ variation becomes even less pronounced after filtering. Due to the smoothing of the CO₂ variation at CDM 6ii, atmospheric CO₂ and δD composition in EDC appear decoupled. This observation seems confirmed when considering the relationship between atmospheric CO₂ change and the duration of NA stadials calculated using tree pollen and the δ¹⁸O composition of planktic foraminifera in an Iberian Margin core (Margari et al., 2010) for MIS 6 (Fig. 3 and Table S4 in the Supplement), or using isotopic records from Greenland ice cores (Rasmussen et al., 2014) for MIS 3 as shown in Fig. 7. Here, the magnitude of atmospheric CO₂ change is generally correlated with the NA stadial duration (r=0.87, n=6) during the early MIS 6 period, and CDM 6ii is by far the shortest of all six events detected in MIS6.

Figure 7 The relationship between NA stadial duration and magnitude of CO₂ increase. Green dots indicate non-Heinrich CDM events during MIS 3, green circles indicate classic Heinrich CDM events during MIS 3 and red-filled green circles indicate non-classic Heinrich CDM events during the MIS 3 period. Blue dots indicate CDM events during MIS 6. Non-classic Heinrich events are defined as ice discharge events with different IRD source signatures from Heinrich events, occurring at transitions in NH ice volume.

We note a similar correlation between the NA stadial duration and atmospheric CO₂ change during MIS 3 (r=0.85, n=14). When the NA stadial duration was shorter than 1500 years, atmospheric CO₂ varied less than 5 ppm (Ahn and Brook, 2014; Bereiter et al., 2012) as is the case for CDM 6ii. Margari et al. (2010) note one exception, AIM 14, during which ice discharge may have led to a stronger perturbation to AMOC.

Both Bereiter et al. (2012) and Ahn and Brook (2014) observed that during short NA stadials which last less than 1300 years, the CO₂ maxima do not appear to have a consistent phase relationship with AIM and CO₂ and δD anomalies are not correlated (Ahn and Brook, 2014). We observe the same to be the case when MIS 6 is included (r=0.07, n=11), however, we cannot exclude that the firn column-induced smoothing at EDC obliterated small atmospheric responses with magnitudes smaller than 5 ppm for such short stadials. In both MIS 3 and 6, CO₂ is highly correlated with δD anomalies during the longer stadials (r=0.78, n=9) with a clear increase in CO₂ (Ahn and Brook, 2014; Bereiter et al., 2012).

We observe that during the last two glacial periods, the amplitude of CO₂ is strongly related to the NA stadial duration (r=0.88, n=20). This observation implies that despite the different climate boundary conditions during the last two glacials, the behavior of atmospheric CO₂ was similar (see Fig. 1) and the overall bipolar seesaw coupling of climate and atmospheric CO₂ acted the same way.

3.4 Leads and lags between CO₂ and the abrupt warming in NH

To better understand the link between the bipolar seesaw mechanism and atmospheric CO₂ variability on millennial timescales, we calculated the varying time lag for each CDM following abrupt warming events in the NH (see Figs. S8–S11 in the Supplement) (Bereiter et al., 2012). Due to the lack of an MIS 6 temperature proxy in Greenland, and due to the difficulty of placing marine temperature proxies (Shackleton et al., 2000) on a precise common chronology with the EDC ice core, in this work CH₄ measurements performed on the EDC ice core were used as a time marker of rapid warming in the NH (Baumgartner et al., 2014; Brook et al., 1996; Huber et al., 2006). Because CH₄ and CO₂ signals are both imprinted in the air bubbles, there is no chronological uncertainty when comparing the timing of changes of those two signals. The only remaining uncertainty is related to analytical uncertainties and to the temporal resolution of the two records. We pick intervals when CH₄ increases rapidly by at least 50 ppb over a time period of less than 1 kyr that correspond to AIM (Buizert et al., 2015; Loulergue et al., 2008). The timing of abrupt CH₄ increases was defined as the midpoint between the beginning of the increase of CH₄ and its peak (see Supplement for details). The age uncertainty of the midpoint is defined by half the time difference between the two endpoints.

Figure 8 shows the shifts of CDM with respect to the onset of the abrupt warming in the NH. During the MIS 6 period, three abrupt NH warmings (as inferred from the CH₄ signal) at 181.5±0.3, 175.4±0.40 and 171.1±0.2 ka (2σ) were found. These events correspond to CDM 6i, CDM 6iii and CDM 6iv, respectively. CDM 6vi, CDM 6v and CDM 6ii do not have corresponding rapid changes in the methane record that fulfill our detection criterion of a 50 ppb increase; this may be due to slow gas trapping as compared to interglacial periods, which could smooth out smaller changes. A synthetic Greenland temperature record (Barker et al., 2011) shows abrupt temperature jumps at CDM 6vi, CDM 6v and CDM 6ii as well. However, this record is calculated using EDC δD, and the large relative chronological uncertainty (∼900 years on average) between ice and gas phases does not allow us to make any conclusions about leads and lags using this record. We therefore exclude these events from our lag analysis.

Figure 8 CDM lags relative to abrupt temperature increases in the NH. Dotted grey lines indicate when climate changes abruptly in the NH as indicated by the CH₄ jumps. During the last glacial period, the AIM number corresponds to the DO number for corresponding DO and AIM events, and for MIS6 numbers correspond to event numbers. (a) Atmospheric CO₂ as measured using the TALDICE ice core during MIS 3, (b) atmospheric CO₂ as measured using the Byrd ice core during MIS 5, (c) atmospheric CO₂ as measured using EDML ice core during MIS 5, (d) new CO₂ composite as derived in this study for MIS 6. For MIS 6, we selected the 3 CDM that correspond to an abrupt methane increase; the other CDM do not correspond to an abrupt change that fulfill our 50 ppb detection criterion (see main text for details). Note that the scale of the y axis is not the same for the four panels.

From MIS 6i to MIS 6iv, the lag of CO₂ with respect to abrupt warmings in the NH, which were identified from this chronological comparison between EDC CH₄ and CO₂, becomes significantly larger. During the earliest MIS 6, atmospheric CO₂ increases rapidly (by ∼4.2 ppm in 240±320 years) following the abrupt CH₄ increase at 181.5±0.3 ka. The peak of CDM 6i is nearly synchronous with the onset of the NH abrupt warming (nonsignificant lag of 240±320 years, Fig. 8). During CDM 6iv and 6iii, CO₂ concentrations show a much slower increase over a duration of ∼3.3 kyr. Here, the CO₂ maximum lags significantly behind the onset of the NH abrupt warming by 1460±270 years and 1110±460 years, respectively (1290±540 years on average, with the error calculated by propagation of the uncertainties in the individual events). Interestingly, these two CDM events occurred during MIS 6d (Fig. 1), when iceberg discharge was muted and the ITCZ is thought to have shifted northward, intensifying monsoon systems in low-latitude Northern Hemisphere regions, such as in Asia, the Apennine Peninsula and the Levant (Ayalon et al., 2002; Bard et al., 2002; Cheng et al., 2016). This may have led to a weaker overturning circulation due to the reduction of the density of the NA surface water, making the AMOC cell shallower with a smaller threshold in NA during MIS 6 than during MIS 3 (Margari et al., 2010). Therefore, the two different CO₂ lag timescales with respect to abrupt warming in NH during MIS 6 might be explained by this difference in background climate conditions. Indeed, these features during MIS 6 appear fully compatible with those observed during the last glacial period (Bereiter et al., 2012). To compare CO₂ variations with respect to abrupt warming in the NH during the last two glacial periods, in this study we also re-estimated the abrupt CH₄ rise during the last glacial period with the same methodologies. Figure 8 also shows the CO₂ evolution during the onset of abrupt warming in the NH during MIS 3 and MIS 5 (Ahn and Brook, 2014; Bereiter et al., 2012). Atmospheric CO₂ during MIS 3, as shown in Fig. 8, was reconstructed from the Talos Dome ice core (TALDICE). For MIS 5 it was obtained from Byrd and the EPICA Dronning Maud Land (EDML) ice core (Ahn and Brook, 2008; Bereiter et al., 2012). In Bereiter et al. (2012), both TALDICE and EDML records were used during MIS 3 and compared to the onset of abrupt warming in the NH as indicated by the co-occurring CH₄ rise. However, here we only use data from TALDICE, which are more accurate due to the narrower age distribution of the gas trapped in the LID (Bereiter et al., 2012). Using the same method as above, the average value of CDM lag with respect to the abrupt warming in NH was calculated. The average CDM lags with respect to the abrupt warming in the NH for the MIS 3 and 5 periods are 770±180 and 280±240 years. Thus, over the course of the last glaciation, the lag of CO₂ maxima with respect to the abrupt NH warming events significantly increased. We observe the same trend through the millennial events depicted during MIS 6, albeit with different absolute lags.

4 Discussion

4.1 Atmospheric CO₂ variability on millennial timescales

Similar to the AIM amplitude (Capron et al., 2010; EPICA Community Members, 2006), we found that the amplitude of atmospheric CO₂ variations is well correlated to the NA stadial duration during MIS 6 and MIS 3, which implies that the amplitude of CO₂ variations might also be affected by the duration of AMOC disruption during the early MIS 6 period (Margari et al., 2010). This hypothesis is also supported by a recent study using oceanic sediment cores from the Southern Ocean (Gottschalk et al., 2020). The authors report that respired carbon levels in the deep South Atlantic decrease when AMOC is weakened during both glacial periods, and the amount of carbon loss in the deep South Atlantic is highly correlated with the duration of NA stadials.

As mentioned above, atmospheric CO₂ on millennial timescales can be controlled by CO₂ exchange between the ocean and the atmosphere, as well as changes of terrestrial carbon stocks. Coupled climate carbon cycle models reported that the variations of atmospheric CO₂ concentration on millennial timescales are mainly dominated by the deep ocean inventory, requiring a few millennia to re-equilibrate to climate change (Schmittner and Galbraith, 2008). On the other hand, the response of the terrestrial biosphere is usually fast (decadal to centennial timescale) (Bouttes et al., 2012; Menviel et al., 2014; Schmittner and Galbraith, 2008). Although different models differ significantly in the CO₂ response to AMOC changes, the initial CO₂ evolution of the terrestrial biosphere and deep ocean to AMOC perturbations are opposite in model simulations (Gottschalk et al., 2019). Thus, due to the opposite direction of CO₂ change of ocean and terrestrial reservoirs, atmospheric CO₂ variations might be muted if the NH duration is short (Bouttes et al., 2012; Menviel et al., 2014; Schmittner and Galbraith, 2008), while during long stadials the carbon release from the ocean dominates. There is, on the other hand, evidence that not all of the processes of CO₂ exchange follow these general trends. For example, atmospheric CO₂ might be changed on centennial timescales by carbon exchange between the deep and surface ocean (Rae et al., 2018) or atmospheric CO₂ might be influenced slowly by soil decomposition (Köhler et al., 2005), and it is important to note that modeling studies do not agree in the amplitude or even the direction of the modeled net CO₂ exchange.

The more prominent CO₂ changes during stadials involved with the Heinrich events may be related to a stronger reduction of the North Atlantic Deep Water (NADW) formation during Heinrich events (Henry et al., 2016; Margari et al., 2010), which would cause a stronger upwelling of deep water in the Southern Ocean (Menviel et al., 2008; Schmittner et al., 2007). These events may reduce stratification in the Southern Ocean due to an increase in salinity of the surface waters and a relative freshening of the deep water (Schmittner et al., 2007). As a result, atmospheric CO₂ can be increased due to upwelling and outgassing of CO₂ in the Southern Ocean (Schmittner et al., 2007; Schmitt et al., 2012). The co-occurring upwelling in the SO during AIM for the last termination has been examined (Anderson et al., 2009), but due to the lack of proxy data with precise age scale for upwelling in the Southern Ocean, this hypothesis cannot be confirmed during MIS 6.

According to the results of Margari et al. (2010), the six AIM events of MIS 6 were likely affected by AMOC perturbations of similar strength. During these events, there is no clear evidence for freshwater perturbation in the NA (Fig. 1), and the strength of the associated AMOC perturbations is estimated to be similar to that during non Heinrich stadials or non-classic Heinrich (different ice-rafted detritus source signatures from Heinrich events, occurring at transitions in ice volume) AIM events of MIS 3. The durations of NA stadials during early MIS 6 (except for AIM 6ii) appear to be longer than those during non-Heinrich AIM events of MIS 3, which might be caused by the different climate boundary conditions during MIS 3 and the increase of the hydrological cycle during the early MIS 6. A longer duration of the AMOC disruption apparently impacts the amplitude of CO₂ variations (Bouttes et al., 2012; Menviel et al., 2008). Considering that the duration of the AMOC disruption may be related with climate background conditions (Bouttes et al., 2012; Menviel et al., 2008), this observation suggests how different climate background conditions may impact atmospheric CO₂. The strength of AMOC perturbations also appears to be an important factor in determining the amplitude of CO₂ variations. For example, the duration of the Heinrich events in MIS 3 (AIM 8, 12 and 14) is shorter than any of the MIS 6 events except for 6ii, but atmospheric CO₂ varied significantly in all three.

The relationship between the amplitude of atmospheric CO₂ variations and the NA stadial duration is explained by the duration of AMOC disruption during the early MIS 6 period. However, the temporal resolution of δ¹⁸O composition of planktonic foraminifera in MD01-2444 and the precision of the age scale are too low to precisely define the duration of stadials during MIS 6. Additional proxy data providing information about climate change in the NH are needed to confirm the relationship between atmospheric CO₂ variations and the NA stadial duration (Fig. 3). The limited available proxy data permit only to formulate a hypothesis for the mechanisms responsible for CO₂ variability during MIS 6 but not to rigorously test it. To compare the behavior of the bipolar seesaw with atmospheric CO₂ variations, additional investigations about AMOC disturbances and their associated climate responses are needed.

4.2 Why did CO₂ lag the abrupt warming in the NH during MIS 6d?

Two different lags of the CO₂ maxima with respect to NH warming are present in the MIS 6 period (Fig. 8). CDM 6i is nearly synchronous with the abrupt warming in the NH (no significant lag of 240±320 years), while the lags for CDM 6iii (1110±460 years) and CDM 6iv (1460±270 years) are much longer. Two modes of CO₂ variations are also observed during the last glacial period. As the last glaciation progressed from MIS 5 to MIS 3 (Figs. 1 and 8), the lag of CO₂ maxima with respect to NH millennial-scale warming significantly increased. This observation may be explained by the different AMOC settings in MIS 5 and MIS 3 (Bereiter et al., 2012). We speculate that, as observed during the last glacial period, the configuration of oceanic circulation during MIS 6d might also be the cause of the change in the time lags between NH abrupt warming events and CO₂ variations during the early MIS 6.

Bereiter et al. (2012) explained that during MIS 3 the oceanic circulation in the Atlantic was in a more “glacial” state, with shallower NADW and carbon-rich Antarctic Bottom Water (AABW) extended to the north, while during MIS 5 ocean circulation was similar to the present, in what can be referred to as a “modern-like” state. At the onset of abrupt warming in the NH, AMOC is thought to accelerate rapidly, delivering heat to the north and resuming the formation of NADW. When the NADW cell expands, AABW is withdrawn and the upwelling of carbon-rich deep water in the Southern Ocean is enhanced. Essentially, over the timescale of ocean overturn, part of the previously expanded carbon-rich southern-sourced water is converted to carbon-poor northern-sourced water.

Thus, during MIS 3 CO₂ continues to be released into the atmosphere for another 500 to 1000 years after the NADW resumption (Bereiter et al., 2012). However, during MIS 5 where AABW was not expanded in the NA, less CO₂ can be released from the ocean to the atmosphere after NADW resumption.

The temporal resolution of proxy data related to oceanic circulation during MIS 3 and 5 is unfortunately not sufficient to validate from the marine realm itself whether the two different modes of CO₂ variations reflect the hypothesized mechanism described above. Modeling studies of the carbon stock in AABW and NADW during MIS 5 and 3 have been attempted. However, dependent on the chosen model, the modes of atmospheric CO₂ variation are different (Gottschalk et al., 2019). Some studies (for example, Menviel et al., 2008) mention an atmospheric CO₂ decrease when the AMOC is resumed. Others, for example Bouttes et al. (2012), confirm an atmospheric CO₂ release from the ocean when AMOC resumed.

In spite of the inconclusive modeling studies (Gottschalk et al., 2019), limited proxy evidence does not exclude the possibility that the configuration of AMOC and its changes over MIS 6 may explain the presence of two different CDM lags. We find this hypothesis to be worth at least a speculative discussion. According to the benthic δ¹³C record in the MD01-2444 core (Margari et al., 2010), the value of benthic δ¹³C during 180–168 ka was lower than during MIS3, which indicates that the NA overturning cell during MIS 6 was likely even shallower than that during MIS 3 (Margari et al., 2010). This implies southern-sourced water masses were more expanded to the north, and the density difference between the northern-sourced water masses and southern-sourced water masses increased. This shallower oceanic circulation during MIS 6 (Margari et al., 2010) could have caused the millennial-scale delays with respect to abrupt NH warming events. It could also explain the longer lag between the abrupt warming in NH and CDM during MIS 6d (1290±540 years on average) when compared to the lags of CDM (770±180 years on average) during MIS 3. However, because of the low accumulation at EDC and its wider age distribution, the estimation of the exact timing of CDM from the EDC ice core might be less accurate compared to that from the TALDICE ice core with a narrower gas age distribution (Bereiter et al., 2012). The remaining uncertainty is related to analytical uncertainties and to the temporal resolution of the two records. To further investigate the exact relationship between CDM and abrupt warming in the NH, additional CO₂ measurements from a higher accumulation site would be helpful.

5 Conclusions

Using new and existing CO₂ data from the EPICA Dome C ice, we reconstruct a high temporal resolution record of atmospheric CO₂ during the MIS 6 period (189–135 ka). In this study, we investigate how different climate background conditions during the last two glacial periods may have impacted atmospheric CO₂. Millennial-scale atmospheric CO₂ changes are revealed during the last two glacial periods, with amplitudes ranging between 15 to 25 ppm, mimicking similar patterns in Antarctic δD variations (Ahn and Brook, 2014; Bereiter et al., 2012). On the other hand, during short NA stadials which last less than 1300 years, atmospheric CO₂ variations are negligible and decoupled from δD in EDC. This finding suggests that during the last two glacial periods the amplitude of millennial CO₂ variations is strongly influenced by the NA stadial duration (r=0.88, n=20).

In the earliest MIS 6 (MIS 6i and 6iv, corresponding to 189 to 169 ka), a change of CO₂ lags with respect to the onset of rapid NH – warming as deduced from atmospheric CH₄ changes – is revealed. CDM 6i (at ∼182 ka) is nearly synchronous with the abrupt warming in the NH (nonsignificant lag of 240±320 years), while the lags during MIS 6d corresponding to CDM 6iv and 6iii (at ∼171 and ∼175 ka, respectively) are much longer, 1290±540 years on average. Similar observations are drawn for the time period covered by our study as for previous studies on MIS 3 and MIS 5 periods, although the lag of CO₂ with respect to NH warming reaches even larger values during MIS 6d. We tentatively attribute this to a generally weaker and shallower AMOC during MIS 6 compared to MIS 3 as suggested by the results from Margari et al. (2010). However, the limited available proxy data from the marine realm only permits an exploratory discussion of the mechanisms responsible for CO₂ variability during MIS 6. Because the boundary conditions of the last glacial period cannot be applied to MIS 6, additional proxy data and multiple modeling studies conducted during the MIS 6 period are needed.