Articles | Volume 16, issue 6
Clim. Past, 16, 2203–2219, 2020
Clim. Past, 16, 2203–2219, 2020

Research article 14 Nov 2020

Research article | 14 Nov 2020

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

Millennial-scale atmospheric CO2 variations during the Marine Isotope Stage 6 period (190–135 ka)
Jinhwa Shin1,a, Christoph Nehrbass-Ahles2,3, Roberto Grilli1, Jai Chowdhry Beeman1, Frédéric Parrenin1, Grégory Teste1, Amaelle Landais4, Loïc Schmidely2, Lucas Silva2, Jochen Schmitt2, Bernhard Bereiter2,5,b, Thomas F. Stocker2, Hubertus Fischer2, and Jérôme Chappellaz1 Jinhwa Shin et al.
  • 1CNRS, Univ. Grenoble-Alpes, Institut des Géosciences de l'Environnement (IGE), Grenoble, France
  • 2Climate and Environmental Physics, Physics Institute, and Oeschger Centre for Climate Change Research, University of Bern, Bern, Switzerland
  • 3Department of Earth Sciences, University of Cambridge, Cambridge, UK
  • 4Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, Gif-sur-Yvette, France
  • 5Laboratory for Air Pollution/Environmental Technology, Empa, Dübendorf, Switzerland
  • acurrent address: Department of Earth and Atmospheric Sciences, University of Alberta, Edmonton, AB, T6G 2E3, Canada
  • bcurrent address: Bruker BioSpin AG, Fällanden, Switzerland

Correspondence: Jérôme Chappellaz (


Using new and previously published CO2 data from the EPICA Dome C ice core (EDC), we reconstruct a new high-resolution record of atmospheric CO2 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 CO2 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 CO2 variation is small (∼5 ppm) and the relationship between temperature variations in EDC and atmospheric CO2 is unclear. The magnitude of CO2 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 CO2 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 CH4 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 CO2 from gas preserved in Antarctic ice sheets (Lüthi et al., 2008; Petit et al., 1999). Comparing atmospheric CO2 records from Antarctic ice cores with proxies of paleoclimate helps us to understand how atmospheric CO2 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 CO2 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 CO2 increased continuously and in parallel to Antarctic temperature increase. Again, at the onset of warming in Greenland, atmospheric CO2 started to decrease (Ahn and Brook, 2008; Bereiter et al., 2012), generally in line with a co-occurring, slow Antarctic cooling. However, the CO2 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 CO2 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 CO2 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 CO2 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 1Proxy 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 δ18Ocalcite 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 CH4 in EDC (dark green dots) (Loulergue et al., 2008) and a composite atmospheric CH4 record from EDC derived in this study (light yellow dots). (f) Composite CO2 record from the EDC ice core derived in this study (orange dots) and a published composite CO2 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 CO2 and climate, we measured 177 new data points of atmospheric CO2 concentrations from the EPICA Dome C (EDC) ice core (7506 S, 12324 E) spanning MIS 6 (190 to 135 ka). This new CO2 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 CO2 dataset measured in 2003. We use all of the available data from the EDC ice core to compile a composite dataset of atmospheric CO2 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 δ15N-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 CO2 concentrations during the early MIS 6, thus establishing a new chronology using new δ15N 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 CH4 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 CO2 data, CH4 measurements are used as a time marker of rapid warming in the NH, as over the last glacial period CH4 and Greenland temperatures are assumed to be synchronous, with a lag of CH4 of less than ∼50 years on average (Baumgartner et al., 2014; Rosen et al., 2014).

2 Methods

2.1CO2 measurements

In this study, we measured atmospheric CO2 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 CO2 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 CO2 measurements were referenced to a secondary gas standard (synthetic air from Air Liquide, Alphagaz 28416000) containing 233.7±0.4 ppm of CO2 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 CO2 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 CO2 should be corrected for contamination caused by the analytical procedure by comparing measured CO2 of the blank tests with the standard gas value. However, it is not feasible to correct CO2 concentrations directly. The CO2 mole fraction is calculated as the ratio between partial pressure of CO2 and total pressure in the measurement line, which is a relative value. Thus, CO2 concentration and the concentration of CO2 contamination are dependent on the total pressure. To avoid this dependence, we corrected the absolute value (partial pressure) of CO2 in the air by the expected partial pressure of CO2 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 CO2 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 CO2 measurements, via four different primary dry air standards for atmospheric CO2, 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 CO2 record was corrected for gravitational fractionation, using the δ15N isotope ratio (Craig et al., 1988). To this end, 88 new data points together with existing δ15 (Landais et al., 2013) covering the late MIS 6 (156.4–139.2 ka) were used. δ15N data were linearly interpolated in age to each corresponding CO2 data point. On average, the correction corresponds to removing 1.2±0.1 ppm from the measured CO2. In total, 177 individual ice samples were measured for CO2 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.2CH4 measurements

We measured the atmospheric CH4 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 CH4 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 CH4 dataset (Loulergue et al., 2008) from EDC has been produced at both IGE and CEP. CH4 measured at CEP used to be systematically higher than CH4 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 δ15N

The water isotopic signature (δD), is, unlike CO2, 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 δ15N 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 δ15N of N2 measurements (Craig et al., 1988; Dreyfus et al., 2010; Sowers et al., 1989) using the following equation:

(1) h = h conv + ( δ 15 N - Ω ( T ) Δ T diff ) Δ m g 1000 R T - 1 .

In this equation, ΔTdiff 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 14N and 15N (kg mol−1), g is the gravitational acceleration (9825 m s−2 for Antarctica) (Parrenin et al., 2013), and R is the universal gas constant (8314 J mol−1K−1). Finally, hconv 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 hconv 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 δ15N 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 Δdepth/z, deduced from the raw δ15N 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 δ15N 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 δ15N 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) Δdepth/z as a function of depth. Red dots represent raw δ15N measurements, and blue dots represent a weighted three-point running mean (see text for details). The dashed vertical grey line indicates Δdepth/z=1. (b) Δdepth (bold line) for the EDC ice core from 1787.5 to 1870.2 m below the surface, deduced from δ15N 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 δ15N 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 δ15N. 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 CO2 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 δ18O 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 δ18O 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 δ18O 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 δ18O 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 3The durations of the six NA stadials during MIS 6: (a) δD of the EDC ice core (Jouzel et al., 2007), (b) synthetic Greenland δ18Oice record (Barker et al., 2011), (c) tree pollen percentage in the MD01-2444 core (Margari et al., 2010) and (d) δ18O 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 δ18O 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 δ18O 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 δ18Oice 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 δ18Oice 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 δ18O 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 δ18Oice 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 CO2 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 CO2 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 CO2 measurements made over the MIS 6 period on the EDC ice core to our new data. First, we compared the two existing CO2 datasets and the new CO2 dataset from EDC (Fig. 4 and Table S5 in the Supplement). There are two published CO2 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 δ15N data in our study. In total, the datasets contain 237 CO2 measurement points. Two samples were excluded because of system operator error. Figure 4 shows CO2 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 CO2 concentrations measured using the ball mill system (Table S5 in the Supplement).

Figure 4Atmospheric CO2 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 CO2 from EDC as measured by the ball mill system (this study). The error bars of the new CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2 from the Vostok ice core (Petit et al., 1999). The error bars of the CO2 data from Vostok (Petit et al., 1999) show the estimated overall accuracy for CO2 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 CO2 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 CO2 record from EDC is shown in Fig. 5.

Figure 5Composite atmospheric CO2 (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 CO2 record for MIS 6

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

Atmospheric CO2 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. CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2 variability, which is the focus of this study.

3.3 Atmospheric CO2 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 CH4 (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, CO2 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 CO2 increased slowly, while during the late MIS 6 period CO2 variation is subdued (Fig. 5).

Figure 6Comparison of climate proxies with atmospheric CO2 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 CH4 in the EDC ice core (Loulergue et al., 2008; this study), (d) δ18O of planktonic foraminifera in the Iberian Margin marine core MD01-2444 (Margari et al., 2010), (e) δ18O 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 CO2 record during the MIS 6 period. The numbers of CDM events are indicated at the top.


While there was an indication of millennial CO2 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 CO2 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 CO2 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 CO2 variations were detected in atmospheric CO2 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 CO2 peak at around 150 ka, representing another potential candidate for a CDM (Fig. 5). This atmospheric CO2 variation is of triangular shape and follows the δD pattern. The change of direction is also associated with a CH4 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 CO2 variation becomes even less pronounced after filtering. Due to the smoothing of the CO2 variation at CDM 6ii, atmospheric CO2 and δD composition in EDC appear decoupled. This observation seems confirmed when considering the relationship between atmospheric CO2 change and the duration of NA stadials calculated using tree pollen and the δ18O 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 CO2 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 7The relationship between NA stadial duration and magnitude of CO2 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 CO2 change during MIS 3 (r=0.85, n=14). When the NA stadial duration was shorter than 1500 years, atmospheric CO2 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 CO2 maxima do not appear to have a consistent phase relationship with AIM and CO2 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, CO2 is highly correlated with δD anomalies during the longer stadials (r=0.78, n=9) with a clear increase in CO2 (Ahn and Brook, 2014; Bereiter et al., 2012).

We observe that during the last two glacial periods, the amplitude of CO2 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 CO2 was similar (see Fig. 1) and the overall bipolar seesaw coupling of climate and atmospheric CO2 acted the same way.

3.4 Leads and lags between CO2 and the abrupt warming in NH

To better understand the link between the bipolar seesaw mechanism and atmospheric CO2 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 CH4 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 CH4 and CO2 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 CH4 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 CH4 increases was defined as the midpoint between the beginning of the increase of CH4 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 CH4 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 8CDM lags relative to abrupt temperature increases in the NH. Dotted grey lines indicate when climate changes abruptly in the NH as indicated by the CH4 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 CO2 as measured using the TALDICE ice core during MIS 3, (b) atmospheric CO2 as measured using the Byrd ice core during MIS 5, (c) atmospheric CO2 as measured using EDML ice core during MIS 5, (d) new CO2 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 CO2 with respect to abrupt warmings in the NH, which were identified from this chronological comparison between EDC CH4 and CO2, becomes significantly larger. During the earliest MIS 6, atmospheric CO2 increases rapidly (by 4.2 ppm in 240±320 years) following the abrupt CH4 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, CO2 concentrations show a much slower increase over a duration of ∼3.3 kyr. Here, the CO2 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 CO2 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 CO2 variations with respect to abrupt warming in the NH during the last two glacial periods, in this study we also re-estimated the abrupt CH4 rise during the last glacial period with the same methodologies. Figure 8 also shows the CO2 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 CO2 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 CH4 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 CO2 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 CO2 variability on millennial timescales

Similar to the AIM amplitude (Capron et al., 2010; EPICA Community Members, 2006), we found that the amplitude of atmospheric CO2 variations is well correlated to the NA stadial duration during MIS 6 and MIS 3, which implies that the amplitude of CO2 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 CO2 on millennial timescales can be controlled by CO2 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 CO2 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 CO2 response to AMOC changes, the initial CO2 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 CO2 change of ocean and terrestrial reservoirs, atmospheric CO2 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 CO2 exchange follow these general trends. For example, atmospheric CO2 might be changed on centennial timescales by carbon exchange between the deep and surface ocean (Rae et al., 2018) or atmospheric CO2 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 CO2 exchange.

The more prominent CO2 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 CO2 can be increased due to upwelling and outgassing of CO2 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 CO2 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 CO2. The strength of AMOC perturbations also appears to be an important factor in determining the amplitude of CO2 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 CO2 varied significantly in all three.

The relationship between the amplitude of atmospheric CO2 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 δ18O 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 CO2 variations and the NA stadial duration (Fig. 3). The limited available proxy data permit only to formulate a hypothesis for the mechanisms responsible for CO2 variability during MIS 6 but not to rigorously test it. To compare the behavior of the bipolar seesaw with atmospheric CO2 variations, additional investigations about AMOC disturbances and their associated climate responses are needed.

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

Two different lags of the CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2 variation are different (Gottschalk et al., 2019). Some studies (for example, Menviel et al., 2008) mention an atmospheric CO2 decrease when the AMOC is resumed. Others, for example Bouttes et al. (2012), confirm an atmospheric CO2 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 δ13C record in the MD01-2444 core (Margari et al., 2010), the value of benthic δ13C 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 CO2 measurements from a higher accumulation site would be helpful.

5 Conclusions

Using new and existing CO2 data from the EPICA Dome C ice, we reconstruct a high temporal resolution record of atmospheric CO2 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 CO2. Millennial-scale atmospheric CO2 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 CO2 variations are negligible and decoupled from δD in EDC. This finding suggests that during the last two glacial periods the amplitude of millennial CO2 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 CO2 lags with respect to the onset of rapid NH – warming as deduced from atmospheric CH4 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 CO2 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 CO2 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.

Data availability

Data available in the Supplement. All data will be available on the NOAA (National Oceanic and Atmospheric Administration) and PANGAEA (Paleoclimatology database websites) shortly.


The supplement related to this article is available online at:

Author contributions

The research was designed by JS, RG, FP and JC. The CO2 measurements were performed by JS with contributions from GT, LS and BB. The data analyses were led by JS with contributions from JCB, RG, FP and JC and HF. The methane data was provided by GT and JC. The nitrogen isotopes data was provided by AL. JS wrote the manuscript with inputs from all authors.

Competing interests

The authors declare that they have no conflict of interest.


This work is a contribution to the “European Project for Ice Coring in Antarctica” (EPICA), a joint European Science Foundation/European Commission scientific program, funded by the European Union and by national contributions from Belgium, Denmark, France, Germany, Italy, the Netherlands, Norway, Sweden, Switzerland and the United Kingdom. The main logistical support was provided by IPEV and PNRA. This is EPICA publication no. 316. It also received funding from the European Community's Seventh Framework Programmes ERC-2011-AdG under grant agreement no. 291062 (ERC ICE&LASERS). As it is part of the PhD work of Jinhwa Shin, it was also supported by the LabEX OSUGat2020 project of the Grenoble Observatory of Sciences of the Universe (OSUG). The Swiss authors also acknowledge long-term financial support for ice core research at the University of Bern by the Swiss National Science Foundation under grants 200020_159563, 200020_172745, 200020_172506 and 20FI21_189533. The authors would like to thank Grégoire Aufresne for providing assistance with the additional methane measurements, Xavier Faïn and Kévin Fourteau for their help with the CH4 analytical system, and Dominique Raynaud for the valuable discussions. We thank Eric Monnin and Urs Siegenthaler for providing additional CO2 data. We would like to thank Julia Gottschalk for the discussions about CO2 variability and climate change during the last two glacial periods.

Financial support

This research has been supported by the ERC-2011-AdG (grant no. 291062, ERC ICE&LASERS), the Swiss National Science Foundation (grant no. 200020_159563), the Swiss National Science Foundation (grant no. 200020_172506), the Grenoble Observatory of Sciences of the Universe (LabEX OSUGat2020 project grant), the Swiss National Science Foundation (grant no. 200020_172745), and the Swiss National Science Foundation (grant no. 20FI21_189533).

Review statement

This paper was edited by Denis-Didier Rousseau and reviewed by four anonymous referees.


Ahn, J. and Brook, E. J.: Atmospheric CO2 and climate on millennial time scales during the last glacial period, Science, 322, 83–85, 2008. 

Ahn, J. and Brook, E. J.: Siple Dome ice reveals two modes of millennial CO2 change during the last ice age, Nat. Commun., 5, 3723,, 2014. 

Anderson, R., Ali, S., Bradtmiller, L., Nielsen, S., Fleisher, M., Anderson, B., and Burckle, L.: Wind-driven upwelling in the Southern Ocean and the deglacial rise in atmospheric CO2, Science, 323, 1443–1448, 2009. 

Ayalon, A., Bar-Matthews, M., and Kaufman, A.: Climatic conditions during marine oxygen isotope stage 6 in the eastern Mediterranean region from the isotopic composition of speleothems of Soreq Cave, Israel, Geology, 30, 303–306, 2002. 

Bard, E., Antonioli, F., and Silenzi, S.: Sea-level during the penultimate interglacial period based on a submerged stalagmite from Argentarola Cave (Italy), Earth Planet. Sci. Lett., 196, 135–146, 2002. 

Barker, S., Knorr, G., Edwards, R. L., Parrenin, F., Putnam, A. E., Skinner, L. C., Wolff, E., and Ziegler, M.: 800 000 years of abrupt climate variability, Science, 334, 347–351, 2011. 

Baumgartner, M., Kindler, P., Eicher, O., Floch, G., Schilt, A., Schwander, J., Spahni, R., Capron, E., Chappellaz, J., Leuenberger, M., Fischer, H., and Stocker, T. F.: NGRIP CH4 concentration from 120 to 10 kyr before present and its relation to a δ15N temperature reconstruction from the same ice core, Clim. Past, 10, 903–920,, 2014. 

Bazin, L., Landais, A., Lemieux-Dudon, B., Toyé Mahamadou Kele, H., Veres, D., Parrenin, F., Martinerie, P., Ritz, C., Capron, E., Lipenkov, V., Loutre, M.-F., Raynaud, D., Vinther, B., Svensson, A., Rasmussen, S. O., Severi, M., Blunier, T., Leuenberger, M., Fischer, H., Masson-Delmotte, V., Chappellaz, J., and Wolff, E.: An optimized multi-proxy, multi-site Antarctic ice and gas orbital chronology (AICC2012): 120–800 ka, Clim. Past, 9, 1715–1731,, 2013. 

Bender, M. L.: Orbital tuning chronology for the Vostok climate record supported by trapped gas composition, Earth Planet. Sc. Lett., 204, 275–289, 2002. 

Bereiter, B., Eggleston, S., Schmitt, J., Nehrbass-Ahles, C., Stocker, T. F., Fischer, H., Kipfstuhl, S., and Chappellaz, J.: Revision of the EPICA Dome C CO2 record from 800 to 600 kyr before present, Geophys. Res. Lett., 42, 542–549, 2015. 

Bereiter, B., Lüthi, D., Siegrist, M., Schüpbach, S., Stocker, T. F., and Fischer, H.: Mode change of millennial CO2 variability during the last glacial cycle associated with a bipolar marine carbon seesaw, Proc. Natl. Acad. Sci., 109, 9755–9760, 2012. 

Bereiter, B., Schwander, J., Lüthi, D., and Stocker, T. F.: Change in CO2 concentration and O2∕N2 ratio in ice cores due to molecular diffusion, Geophys. Res. Lett., 36,, 2009. 

Bereiter, B., Stocker, T. F., and Fischer, H.: A centrifugal ice microtome for measurements of atmospheric CO2 on air trapped in polar ice cores, Atmos. Meas. Tech., 6, 251–262,, 2013. 

Berger, A. L.: Long-Term Variations of Caloric Insolation Resulting from the Earth's Orbital Elements 1, Quat. Res., 9, 139–167, 1978. 

Blunier, T. and Brook, E. J.: Timing of millennial-scale climate change in Antarctica and Greenland during the last glacial period, Science, 291, 109–112, 2001. 

Bond, G., Heinrich, H., Broecker, W., Labeyrie, L., McManus, J., Andrews, J., Huon, S., Jantschik, R., Clasen, S., and Simet, C.: Evidence for massive discharges of icebergs into the North Atlantic ocean during the last glacial period, Nature, 360, 245–249, 1992. 

Bouttes, N., Roche, D. M., and Paillard, D.: Systematic study of the impact of fresh water fluxes on the glacial carbon cycle, Clim. Past, 8, 589–607,, 2012. 

Bozbiyik, A., Steinacher, M., Joos, F., Stocker, T. F., and Menviel, L.: Fingerprints of changes in the terrestrial carbon cycle in response to large reorganizations in ocean circulation, Clim. Past, 7, 319–338,, 2011. 

Bréant, C., Landais, A., Orsi, A., Martinerie, P., Extier, T., Prié, F., Stenni, B., Jouzel, J., Masson-Delmotte, V., and Leuenberger, M.: Unveiling the anatomy of Termination 3 using water and air isotopes in the Dome C ice core, East Antarctica, Quaternary Sci. Rev., 211, 156–165, 2019. 

Broecker, W., Bond, G., Klas, M., Clark, E., and McManus, J.: Origin of the northern Atlantic's Heinrich events, Clim. Dyn., 6, 265–273, 1992. 

Brook, E. J., Sowers, T., and Orchardo, J.: Rapid variations in atmospheric methane concentration during the past 110 000 yrs, Science, 273, 1087–1091, 1996. 

Buizert, C., Adrian, B., Ahn, J., Albert, M., Alley, R. B., Baggenstos, D., Bauska, T. K., Bay, R. C., Bencivengo, B. B., Bentley, C. R., Brook, E. J., Chellman, N. J., Clow, G. D., Cole-Dai, J., Conway, H., Cravens, E., Cuffey, K. M., Dunbar, N. W., Edwards, J. S., Fegyveresi, J. M., Ferris, D. G., Fitzpatrick, J. J., Fudge, T. J., Gibson, C. J., Gkinis, V., Goetz, J. J., Gregory, S., Hargreaves, G. M., Iverson, N., Johnson, J. A., Jones, T. R., Kalk, M. L., Kippenhan, M. J., Koffman, B. G., Kreutz, K., Kuhl, T. W., Lebar, D. A., Lee, J. E., Marcott, S. A., Markle, B. R., Maselli, O. J., McConnell, J. R., McGwire, K. C., Mitchell, L. E., Mortensen, N. B., Neff, P. D., Nishiizumi, K., Nunn, R. M., Orsi, A. J., Pasteris, D. R., Pedro, J. B., Pettit, E. C., Price, P. B., Priscu, J. C., Rhodes, R. H., Rosen, J. L., Schauer, A. J., Schoenemann, S. W., Sendelbach, P. J., Severinghaus, J. P., Shturmakov, A. J., Sigl, M., Slawny, K. R., Souney, J. M., Sowers, T. A., Spencer, M. K., Steig, E. J., Taylor, K. C., Twickler, M. S., Vaughn, B. H., Voigt, D. E., Waddington, E. D., Welten, K. C., Wendricks, A. W., White, J. W. C., Winstrup, M., Wong, G. J., and Woodruff, T. E.: Precise interpolar phasing of abrupt climate change during the last ice age, Nature, 520, 661–665,, 2015. 

Capron, E., Landais, A., Lemieux-Dudon, B., Schilt, A., Masson-Delmotte, V., Buiron, D., Chappellaz, J., Dahl-Jensen, D.,Johnsen, S., Leuenberger, M., Loulergue, L., and Oerter, H.: Synchronising EDML and North GRIP ice cores using δ18O of atmospheric oxygen (δ18Oatm) and CH4 measurements over MIS5 (80–123 kyr), Quaternary Sci. Rev., 29, 222–234,, 2010. 

Cheng, H., Edwards, R. L., Sinha, A., Spötl, C., Yi, L., Chen, S., Kelly, M., Kathayat, G., Wang, X., and Li, X.: The Asian monsoon over the past 640 000 years and ice age terminations, Nature, 534, 640–646,, 2016. 

Craig, H., Horibe, Y., and Sowers, T.: Gravitational separation of gases and isotopes in polar ice caps, Science, 242, 1675–1678, 1988. 

de Abreu, L., Shackleton, N. J., Schönfeld, J., Hall, M., and Chapman, M.: Millennial-scale oceanic climate variability off the Western Iberian margin during the last two glacial periods, Mar. Geol., 196, 1–20, 2003. 

Denton, G. H., Anderson, R. F., Toggweiler, J., Edwards, R., Schaefer, J., and Putnam, A.: The last glacial termination, Science, 328, 1652–1656, 2010. 

Dreyfus, G. B., Jouzel, J., Bender, M. L., Landais, A., Masson-Delmotte, V., and Leuenberger, M.: Firn processes and δ15N: potential for a gas-phase climate proxy, Quaternary Sci. Rev., 29, 28–42, 2010. 

EPICA Community Members: One-to-one coupling of glacial climate variability in Greenland and Antarctica, Nature, 444, 195–198, 2006. 

Etheridge, D., Pearman, G., and Fraser, P.: Changes in tropospheric methane between 1841 and 1978 from a high accumulation-rate Antarctic ice core, Tellus B, 44, 282–294, 1992. 

Fischer, H., Schmitt, J., Lüthi, D., Stocker, T. F., Tschumi, T.,Parekh, P., Joos, F., Köhler, P., Völker, C., Gersonde, R., Barbante, C., Le Floch, M., Raynaud, D., and Wolff, E.: The role of Southern Ocean processes in orbital and millennial CO2 variations – A synthesis, Quat. Sci. Rev., 29, 193–205, 2010. 

Fourteau, K., Faïn, X., Martinerie, P., Landais, A., Ekaykin, A. A., Lipenkov, V. Ya., and Chappellaz, J.: Analytical constraints on layered gas trapping and smoothing of atmospheric variability in ice under low-accumulation conditions, Clim. Past, 13, 1815–1830,, 2017. 

Gottschalk, J., Battaglia, G., Fischer, H., Frölicher, T. L., Jaccard, S. L., Jeltsch-Thömmes, A., Joos, F., Köhler, P., Meissner, K. J., and Menviel, L., Nehrbass-Ahles, C., Schmitt, J., Schmittner, A., Skinner, L. C., and Stocker, T.: Mechanisms of millennial-scale atmospheric CO2 change in numerical model simulations, Quaternary Sci. Rev.,220, 30–74, 2019. 

Gottschalk, J., Skinner, L. C., Jaccard, S. L., Menviel, L., Nehrbass-Ahles, C., and Waelbroeck, C.: Southern Ocean link between changes in atmospheric CO2 levels and northern-hemisphere climate anomalies during the last two glacial periods, Quaternary Sci. Rev., 230, 106067,, 2020. 

Grachev, A. M. and Severinghaus, J. P.: Laboratory determination of thermal diffusion constants for 29N228N2 in air at temperatures from −60 to 0 C for reconstruction of magnitudes of abrupt climate changes using the ice core fossil–air paleothermometer, Geochim. Cosmochim. Ac., 67, 345–360, 2003. 

Heinrich, H.: Origin and consequences of cyclic ice rafting in the northeast Atlantic Ocean during the past 130 000 years, Quaternary Res., 29, 142–152, 1988. 

Henry, L., McManus, J. F., Curry, W. B., Roberts, N. L., Piotrowski, A. M., and Keigwin, L. D.: North Atlantic ocean circulation and abrupt climate change during the last glaciation, Science, 353, 470–474, 2016. 

Higgins, J. A., Kurbatov, A. V., Spaulding, N. E., Brook, E., Introne, D. S., Chimiak, L. M., Yan, Y., Mayewski, P. A., and Bender, M. L.: Atmospheric composition 1 million years ago from blue ice in the Allan Hills, Antarctica, Proc. Natl. Acad. Sci., 112, 6887–6891, 201420232,, 2015. 

Huber, C., Leuenberger, M., Spahni, R., Flückiger, J., Schwander, J., Stocker, T. F., Johnsen, S., Landais, A., and Jouzel, J.: Isotope calibrated Greenland temperature record over Marine Isotope Stage 3 and its relation to CH4, Earth Planet. Sci. Lett., 243, 504–519, 2006. 

Jouzel, J., Masson-Delmotte, V., Cattani, O., Dreyfus, G., Falourd, S., Hoffmann, G., Minster, B., Nouet, J., Barnola, J. M., Chappellaz, J., Fischer, H., Gallet, J. C., Johnsen, S., Leuenberger, M., Loulergue, L., Luethi, D., Oerter, H., Parrenin, F., Raisbeck, G., Raynaud, D., Schilt, A., Schwander, A., Selmo, E., Souchez, R., Spahni, R., Stauffer, B., Steffensen, J. P., Stenni, B., Stocker, T. F., Tison, J. L., Werner, M., and Wolff, E. W.: Orbital and Millennial Antarctic Climate Variability over the Past 800 000 Years, Science 317, 793–796,, 2007. 

Kawamura, K., Abe-Ouchi, A., Motoyama, H., Ageta, Y., Aoki, S., Azuma, N., Fujii, Y., Fujita, K., Fujita, S., and Fukui, K.: State dependence of climatic instability over the past 720 000 years from Antarctic ice cores and climate modeling, Sci. Adv., 3, e1600446,, 2017. 

Köhler, P., Joos, F., Gerber, S., and Knutti, R.: Simulated changes in vegetation distribution, land carbon storage, and atmospheric CO2 in response to a collapse of the North Atlantic thermohaline circulation, Clim. Dyn., 25, 689,, 2005. 

Krinner, G., Raynaud, D., Doutriaux, C., and Dang, H.: Simulations of the Last Glacial Maximum ice sheet surface climate: Implications for the interpretation of ice core air content, J. Geophys. Res., 105, 2059–2070, 2000. 

Landais, A., Barnola, J., Kawamura, K., Caillon, N., Delmotte, M., Ommen, T. V., Dreyfus, G., Jouzel, J., Masson-Delmotte, V., Minster, B., Freitag, J., Leuenberger, M., Schwander, J., Huber, C., Etheridge, D., and Morgan, V.: Firn-air δ15N in modern polar sites and glacial-interglacial ice: a model-data mismatch during glacial periods in Antarctica?, Quaternary Sci. Rev., 25, 49–62, 2006. 

Landais, A., Dreyfus, G., Capron, E., Jouzel, J., Masson-Delmotte, V., Roche, D. M., Prié, F., Caillon, N., Chappellaz, J., Leuenberger, M., Lourantou, A., Parrenin, F., Raynaud, D., and Teste, G.: Two-phase change in CO2, Antarctic temperature and global climate during Termination II, Nat. Geosci., 6, 1062–1065,, 2013. 

Loulergue, L., Parrenin, F., Blunier, T., Barnola, J.-M., Spahni, R., Schilt, A., Raisbeck, G., and Chappellaz, J.: New constraints on the gas age-ice age difference along the EPICA ice cores, 0–50 kyr, Clim. Past, 3, 527–540,, 2007. 

Loulergue, L., Schilt, A., Spahni, R., Masson-Delmotte, V., Blunier, T., Lemieux, B., Barnola, J.-M., Raynaud, D., Stocker, T. F., and Chappellaz, J.: Orbital and millennial-scale features of atmospheric CH4 over the past 800 000 years, Nature, 453, 383,, 2008.  

Lourantou, A., Chappellaz, J., Barnola, J.-M., Masson-Delmotte, V., and Raynaud, D.: Changes in atmospheric CO2 and its carbon isotopic ratio during the penultimate deglaciation, Quaternary Sci. Rev., 29, 1983–1992, 2010. 

Lüthi, D., Bereiter, B., Stauffer, B., Winkler, R., Schwander, J., Kindler, P., Leuenberger, M., Kipfstuhl, S., Capron, E., and Landais, A.: CO2 and O2∕N2 variations in and just below the bubble-clathrate transformation zone of Antarctic ice cores, Earth Planet. Sci. Lett., 297, 1–2, 226–233,, 2010. 

Luthi, D., Le Floch, M., Bereiter, B., Blunier, T., Barnola, J.-M., ¨ Siegenthaler, U., Raynaud, D., Jouzel, J., Fischer, H., Kawamura, K., and Stocker, T. F.: High-resolution carbon dioxide concentration record 650 000–800 000 years before present, Nature, 453, 379–382,, 2008. 

Marcott, S. A., Bauska, T. K., Buizert, C., Steig, E. J., Rosen, J. L., Cuffey, K. M., Fudge, T. J., Severinghaus, J. P., Ahn, J., Kalk, M. L., McConnell, J. R., Sowers, T., Taylor, K. C., White, J. W., and Brook, E. J.: Centennial scale changes in the global carbon cycle during the last deglaciation, Nature, 514, 616–619,, 2014. 

Margari, V., Skinner, L., Tzedakis, P., Ganopolski, A., Vautravers, M., and Shackleton, N.: The nature of millennial-scale climate variability during the past two glacial periods, Nat. Geosci., 3, 127–131,, 2010. 

McManus, J. F., Anderson, R. F., Broecker, W. S., Fleisher, M. Q., and Higgins, S. M.: Radiometrically determined sedimentary fluxes in the sub-polar North Atlantic during the last 140 000 years, Earth Planet. Sci. Lett., 155, 29–43, 1998. 

McManus, J. F., Oppo, D. W., and Cullen, J. L.: A 0.5 million-year record of millennial-scale climate variability in the North Atlantic, Science, 283, 971–975, 1999. 

Menviel, L., England, M. H., Meissner, K., Mouchet, A., and Yu, J.: Atlantic-Pacific seesaw and its role in outgassing CO2 during Heinrich events, Paleoceanography, 29, 58–70, 2014. 

Menviel, L., Timmermann, A., Mouchet, A., and Timm, O.: Meridional reorganizations of marine and terrestrial productivity during Heinrich events, Paleoceanography, 23, PA1203,, 2008. 

Monnin, E., Steig, E. J., Siegenthaler, U., Kawamura, K., Schwander, J., Stauffer, B., Stocker, T. F., Morse, D. C., Barnola, J.-M., Bellier, B., Raynaud, D., and Fischer, H.: Evidence for substantial accumulation rate variability in Antarctica during the Holocene through synchronization of CO2 in the Taylor Dome, Dome C and DML ice cores, Earth Planet. Sc. Lett., 224, 45–54, 2004 

Nehrbass-Ahles, C., Shin, J., Schmitt, J., Bereiter, B., Joos, F., Schilt, A., Schmidely, L., Silva, L., Teste, G., Grilli, R., Chappellaz, J., Hodell, D., Fischer, H., and Stocker, T. F.: Abrupt CO2 release to the atmosphere under both glacial and early interglacial conditions, Science, 369, 1000–1005, 2020. 

Parrenin, F., Barker, S., Blunier, T., Chappellaz, J., Jouzel, J., Landais, A., Masson-Delmotte, V., Schwander, J., and Veres, D.: On the gas-ice depth difference (δdepth) along the EPICA Dome C ice core, Clim. Past, 8, 1239–1255,, 2012. 

Parrenin, F., Barnola, J.-M., Beer, J., Blunier, T., Castellano, E., Chappellaz, J., Dreyfus, G., Fischer, H., Fujita, S., Jouzel, J., Kawamura, K., Lemieux-Dudon, B., Loulergue, L., Masson-Delmotte, V., Narcisi, B., Petit, J.-R., Raisbeck, G., Raynaud, D., Ruth, U., Schwander, J., Severi, M., Spahni, R., Steffensen, J. P., Svensson, A., Udisti, R., Waelbroeck, C., and Wolff, E.: The EDC3 chronology for the EPICA Dome C ice core, Clim. Past, 3, 485–497,, 2007. 

Parrenin, F., Masson-Delmotte, V., Köhler, P., Raynaud, D., Paillard, D., Schwander, J., Barbante, C., Landais, A., Wegner, A., and Jouzel, J.: Synchronous change of atmospheric CO2 and Antarctic temperature during the last deglacial warming, Science, 339, 1060–1063, 2013. 

Pedro, J. B., Jochum, M., Buizert, C., He, F., Barker, S., and Rasmussen, S. O.: Beyond the bipolar seesaw: Toward a process understanding of interhemispheric coupling, Quaternary Sci. Rev., 192, 27–46, 2018. 

Petit, J. R., Jouzel, J., Raynaud, D., Barkov, N. I., Barnola, J.-M., Basile, I., Bender, M., Chapellaz, J., Davis, M., Delaygue, G., Delmotte, M., Kotlyakov, V. M., Legrand, M., Lipenkov, V. Y., Lorius, C., Pepin, L., Ritz, C., Saltzmann, E., and Stievenard, M.: Climate and atmospheric history of the past 420 000 years from the Vostok ice core, Antarctica, Nature, 399, 429–436,, 1999. 

Rae, J. W. B., Burke, A., Robinson, L. F., Adkins, J. F., Chen, T., Cole, C., Greenop, R., Li, T., Littley, E. F. M., Nita, D. C., and Stewart, J. A.: CO2 storage and release in the deep Southern Ocean on millennial to centennial timescales, Nature, 562, 569–573, 2018. 

Railsback, L. B., Gibbard, P. L., Head, M. J., Voarintsoa, N. R. G., and Toucanne, S.: An optimized scheme of lettered marine isotope substages for the last 1.0 million years, and the climatostratigraphic nature of isotope stages and substages, Quaternary Sci. Rev., 111, 94–106, 2015. 

Rasmussen, S. O., Bigler, M., Blockley, S. P., Blunier, T., Buchardt, S. L., Clausen, H. B., Cvijanovic, I., Dahl-Jensen, D., Johnsen, S. J., Fischer, H., Gkinis, V., Guillevic, M., Hoek, W. Z., Lowe, J. J., Pedro, J. B., Popp, T., Seierstad, I. K., Steffensen, J. P., Svensson, A. M., Vallelonga, P., Vinther, B. M., Walker, M. J. C., Wheatley, J. J., and Winstrup, M.: A stratigraphic framework for abrupt climatic changes during the Last Glacial period based on three synchronized Greenland ice-core records: refining and extending the INTIMATE event stratigraphy, Quaternary Sci. Rev., 106, 14–28, 2014. 

Rosen, J. L., Brook, E. J., Severinghaus, J. P., Blunier, T., Mitchell, L. E., Lee, J. E., Edwards, J. S., and Gkinis, V.: An ice core record of near-synchronous global climate changes at the Bølling transition, Nat. Geosci., 7, 459,, 2014. 

Schaefer, H., Lourantou, A., Chappellaz, J., Lüthi, D., Bereiter, B., and Barnola, J.-M.: On the suitability of partially clathrated ice for analysis of concentration and δ13C of palaeo-atmospheric CO2, Earth Planet. Sci. Lett., 307, 334–340, 2011. 

Schmitt, J., Schneider, R., Elsig, J., Leuenberger, D., Lourantou, A., Chappellaz, J., Köhler, P., Joos, F., Stocker, T. F., and Leuenberger, M.: Carbon isotope constraints on the deglacial CO2 rise from ice cores, Science, 336, 711–714, 2012. 

Schmittner, A., Brook, E. J., and Ahn, J.: Impact of the ocean's overturning circulation on atmospheric CO2, in: Ocean Circulation: Mechanisms and Impacts, AGU Geophysical Monograph Series, 173, American Geophysical Union, Washington DC, 315–334, 2007. 

Schmittner, A. and Galbraith, E. D.: Glacial greenhouse-gas fluctuations controlled by ocean circulation changes, Nature, 456, 373,, 2008. 

Schneider, R., Schmitt, J., Köhler, P., Joos, F., and Fischer, H.: A reconstruction of atmospheric carbon dioxide and its stable carbon isotopic composition from the penultimate glacial maximum to the last glacial inception, Clim. Past, 9, 2507–2523,, 2013. 

Shackleton, N. J., Hall, M. A., and Vincent, E.: Phase relationships between millennial-scale events 64 000–24 000 years ago, Paleoceanography, 15, 565–569, 2000. 

Siegenthaler, U., Monnin, E., Kawamura, K., Spahni, R., Schwander, J., Stauffer, B., Stocker, T. F., Barnola, J.-M., and Fischer, H.: Supporting evidence from the EPICA Dronning Maud Land ice core for atmospheric CO2 changes during the past millennium, Tellus B, 57, 51–57, 2005. 

Sigman, D. M. and Boyle, E. A.: Glacial/interglacial variations in atmospheric carbon dioxide, Nature, 407, 859,, 2000. 

Sigman, D. M., De Boer, A. M., and Haug, G. H.: Antarctic stratification, atmospheric water vapor, and Heinrich events: A hypothesis for late Pleistocene deglaciations, 173, American Geophysical Union, Washington DC, 315–334, 2007. 

Sowers, T., Bender, M., and Raynaud, D.: Elemental and isotopic composition of occluded O2 and N2 in polar ice, J. Geophys. Res., 94, 5137–5150, 1989. 

Spahni, R., Chappellaz, J., Stocker, T., Loulergue, L., Hausmmann, G., Kawamura, K., Flückiger, J., Schwander, J., Raynaud, D., Masson-Delmotte, V., and Jouzel, J.: Atmospheric Methane and Nitrous Oxide of the Late Pleistocene from Antarctic Ice Cores, Science, 310, 1317–1321,, 2005. 

Stocker, T. F. and Johnsen, S. J.: A minimum thermodynamic model for the bipolar seesaw, Paleoceanography, 18, PA000920,, 2003. 

Tans, P. P., Crotwell, A. M., and Thoning, K. W.: Abundances of isotopologues and calibration of CO2 greenhouse gas measurements, Atmos. Meas. Tech., 10, 2669–2685,, 2017. 

Toggweiler, J. R., Russell, J. L., and Carson, S. R.: Midlatitude westerlies, atmospheric CO2, and climate change during the ice ages, Paleoceanography, 21, PA2005,, 2006. 

Short summary
We reconstruct atmospheric CO2 from the EPICA Dome C ice core during Marine Isotope Stage 6 (185–135 ka) to understand carbon mechanisms under the different boundary conditions of the climate system. The amplitude of CO2 is highly determined by the Northern Hemisphere stadial duration. Carbon dioxide maxima show different lags with respect to the corresponding abrupt CH4 jumps, the latter reflecting rapid warming in the Northern Hemisphere.