Articles | Volume 15, issue 3
Research article
22 May 2019
Research article |  | 22 May 2019

Antarctic temperature and CO2: near-synchrony yet variable phasing during the last deglaciation

Jai Chowdhry Beeman, Léa Gest, Frédéric Parrenin, Dominique Raynaud, Tyler J. Fudge, Christo Buizert, and Edward J. Brook

The last deglaciation, which occurred from 18 000 to 11 000 years ago, is the most recent large natural climatic variation of global extent. With accurately dated paleoclimate records, we can investigate the timings of related variables in the climate system during this major transition. Here, we use an accurate relative chronology to compare temperature proxy data and global atmospheric CO2 as recorded in Antarctic ice cores. In addition to five regional records, we compare a δ18O stack, representing Antarctic climate variations with the high-resolution robustly dated WAIS Divide CO2 record (West Antarctic Ice Sheet). We assess the CO2 and Antarctic temperature phase relationship using a stochastic method to accurately identify the probable timings of changes in their trends. Four coherent changes are identified for the two series, and synchrony between CO2 and temperature is within the 95 % uncertainty range for all of the changes except the end of glacial termination 1 (T1). During the onset of the last deglaciation at 18 ka and the deglaciation end at 11.5 ka, Antarctic temperature most likely led CO2 by several centuries (by 570 years, within a range of 127 to 751 years, 68 % probability, at the T1 onset; and by 532 years, within a range of 337 to 629 years, 68 % probability, at the deglaciation end). At 14.4 ka, the onset of the Antarctic Cold Reversal (ACR) period, our results do not show a clear lead or lag (Antarctic temperature leads by 50 years, within a range of −137 to 376 years, 68 % probability). The same is true at the end of the ACR (CO2 leads by 65 years, within a range of 211 to 117 years, 68 % probability). However, the timings of changes in trends for the individual proxy records show variations from the stack, indicating regional differences in the pattern of temperature change, particularly in the WAIS Divide record at the onset of the deglaciation; the Dome Fuji record at the deglaciation end; and the EDML record after 16 ka (EPICA Dronning Maud Land, where EPICA is the European Project for Ice Coring in Antarctica). In addition, two changes – one at 16 ka in the CO2 record and one after the ACR onset in three of the isotopic temperature records – do not have high-probability counterparts in the other record. The likely-variable phasing we identify testify to the complex nature of the mechanisms driving the carbon cycle and Antarctic temperature during the deglaciation.

1 Introduction

Glacial–interglacial transitions, or deglaciations, mark the paleorecord approximately every 100 000 years over the past million years or so (Williams et al.1997; Lisiecki and Raymo2005; Jouzel et al.2007). The last deglaciation, often referred to as glacial termination 1 (T1), offers a case study for a large global climatic change, very likely in the 3–8 C range on the regional scale (Masson-Delmotte et al.2013), and thought to be initiated by an orbitally driven insolation forcing (Hays et al.1976; Berger1978; Kawamura et al.2007). The canonical interpretation of this apparent puzzle is that insolation acts as a pacemaker of climatic cycles and the amplitude of glacial–interglacial transitions is mainly driven by two strong climatic feedbacks: atmospheric CO2 and continental-ice-surface–albedo changes. However, the mechanisms that control the CO2 rise are still a matter of debate. Accordingly, reconstructing the phase relationship (leads and lags) between climate variables and CO2 during the last termination has become important, and has a substantial history in ice core research (Barnola et al.1991; Raynaud and Siegenthaler1993; Caillon et al.2003; Pedro et al.2012; Parrenin et al.2013).

Global temperature has been shown to lag CO2 on average during T1 (Shakun et al.2012), supporting the importance of CO2 as an amplifier of orbitally driven global-scale warming. But Antarctic temperature and CO2 concentrations changed much more coherently as T1 progressed. Indeed, midway through the glacial–interglacial transition, Antarctic warming and CO2 increase slowed and even reversed during a period of about 2000 years, coinciding with a warm period in the north called the Bølling–Allerød (BA). The respective period of cooling in Antarctica is called the Antarctic Cold Reversal (ACR). A period of cooling in the Northern Hemisphere known as the Younger Dryas (YD), followed the BA, coinciding with a period of warming in the Southern Hemisphere (SH).

High-latitude Southern Hemisphere paleotemperature series – including Southern Ocean temperature – varied similarly to Antarctic temperature during T1 (Shakun et al.2012; Pedro et al.2016). Upwelling from the Southern Ocean is thought to have played an important role in the deglacial CO2 increases (Anderson et al.2009; Burke and Robinson2012; Schmitt et al.2012; Rae et al.2018). The Atlantic Meridional Overturning Circulation, or AMOC, a major conduit of heat between the Northern Hemisphere and Southern Hemisphere and component of the bipolar seesaw – the umbrella term encompassing the mechanisms thought to control the seemingly alternating variations in Northern Hemisphere and Southern Hemisphere temperature (Stocker and Johnsen2003; Pedro et al.2018) – is thought to have influenced Southern Ocean upwelling during the deglaciation (Anderson et al.2009; Skinner et al.2010). A weakening of the oceanic biological carbon pump appears to have dominated the deglacial CO2 increase until 15.5 ka (kiloannum before 1950 CE), when rising ocean temperature likely began to play a role as well (Bauska et al.2016).

Ice sheets are exceptional archives of past climates and atmospheric composition. Local temperature is recorded in the isotopic composition of snow/ice (NorthGRIP Project Members2004; Jouzel et al.2007) due to the temperature-dependent fractionation of water isotopes (Lorius and Merlivat1975; Johnsen et al.1989). The concentration of continental dust in ice sheets is a proxy of continental aridity, atmospheric transport intensity and precipitation. Finally, air bubbles enclosed in ice sheets are near-direct samples of the past atmosphere. However, the age of the air bubbles is younger than the age of the surrounding ice, since air is locked-in at the base of the firn (on the order of 70 m below the surface on the West Antarctic Ice Sheet (WAIS) Divide) at the lock-in depth (LID) (Buizert and Severinghaus2016). The firn, from top to bottom, is composed of a convective zone (CZ), where the air is mixed vigorously, and a diffusive zone (DZ), where molecular diffusion dominates transport. Firn densification models can be used to estimate the LID and the corresponding age difference (Sowers et al.1992).

Atmospheric CO2 concentrations, recorded in the air bubbles enclosed in ice sheets, are better preserved in Antarctic ice than in Greenland ice because the latter has much higher concentrations of organic material and carbonate dust (Raynaud et al.1993; Anklin et al.1995). Measured on the Vostok and EPICA Dome C ice cores (EPICA is the European Project for Ice Coring in Antarctica), the long-term history of CO2 (Lüthi et al.2008) covers the last 800 kyr.

Early studies suggested that at the initiation of the termination around 18 ka, just after the Last Glacial Maximum (LGM), Antarctic temperature started to warm 800±600 years before CO2 began to increase (Monnin et al.2001), a result that was sometimes misinterpreted to mean that CO2 was not an important amplification factor of the deglacial temperature increase. This study used measurements from the EPICA Dome C (EDC) ice core (Jouzel et al.2007) and a firn densification model to determine the air chronology. However, this firn densification model was later shown to be in error by several centuries for low-accumulation sites such as EDC during glacial periods (Loulergue et al.2007; Parrenin et al.2012).

Two more recent works (Pedro et al.2012; Parrenin et al.2013) used stacked temperature records and improved estimates of the age difference between ice and air records to more accurately estimate the relative timing of changes in Antarctic temperature and atmospheric CO2 concentration. In the first of these studies, measurements from the higher accumulation ice cores at Siple Dome and Byrd Station were used to decrease the uncertainty in the ice–air age shift, and indicated that CO2 lagged Antarctic temperature by 0–400 years on average during the last deglaciation (Pedro et al.2012). The second study (Parrenin et al.2013) used measurements from the low accumulation EDC ice core but circumvented the use of firn densification models by using the nitrogen isotope ratio δ15N of N2 as a proxy of the DZ height, assuming that the height of the CZ was negligible during the study period. CO2 and Antarctic temperature were found to be in phase at the beginning of TI (-10±160 years) and at the end of the ACR period (-60±120 years), but CO2 was found to lag Antarctic temperature by several centuries at the beginning of the Antarctic Cold Reversal (260±130 years) and at the end of the deglacial warming in Antarctica (500±90 years). The end of the deglacial warming in Antarctica occurred roughly 2 centuries after the onset of the Holocene period dated at 11.7 ka according to the International Commission on Stratigraphy. However, the assumption that δ15N reflects DZ height is imperfect as it may underestimate the DZ height for sites with strong barometric pumping and layering (Buizert and Severinghaus2016).

Figure 1Drilling locations of the ice cores where the CO2 and isotopic paleotemperature records included in this study were measured.


A new CO2 record of unprecedented high resolution (Marcott et al.2014) from the WAIS Divide (WD) ice core merits the reopening of this investigation. The air chronology of WAIS Divide is well constrained thanks to a relatively high accumulation rate and to accurate nitrogen-15 measurements (Buizert et al.2015). The WAIS record evidences centennial-scale changes in the global carbon cycle during the last deglaciation superimposed on more gradual, millennial-scale trends that bear a resemblance to Antarctic temperature (Marcott et al.2014).

The deglacial temperature rise seen at WD is structurally similar to that at other Antarctic sites. However, the rise in West Antarctic temperature shows an early warming starting around 21 ka, following local insolation (Cuffey et al.2016). This early warming trend is much more gradual in records from East Antarctic ice cores. The difference between the two records may be related to sea ice conditions around East Antarctica and West Antarctica, and perhaps to elevation changes (WAIS Divide Project Members2013; Cuffey et al.2016). The temperature record at WAIS Divide shows an acceleration in warming around 18 ka which is also present in East Antarctic records (WAIS Divide Project Members2013).

On the much shorter timescales of the satellite era, Jones et al. (2016) note differing temperature trends at the drilling sites of the five cores used in this study. On the other hand, the interpretation of individual isotopic records can prove to be complicated, as local effects – including those of ice sheet elevation change and sea ice extent – are difficult to correct.

In the present work we refine our knowledge of leads and lags between Antarctic temperature and CO2. We use a stack of accurately synchronized Antarctic temperature records (Buizert et al.2018) to reduce local signals, placed using volcanic matching on the WAIS Divide chronology (WD2014). We then compare the temperature stack to the high-resolution WAIS Divide CO2 record by determining the probable timings of changes in trends, and calculate probable change point timings for the five individual isotope-derived records used in the stack as well.

2 Methods and data

2.1 Temperature stack and ice chronology

We use the δ18O stack developed by Buizert et al. (2018) (referred to hereafter as Antarctic temperature Stack 3, or ATS3) to represent Antarctic temperature. The use of the stack allows us to remove local influences and noise in the individual records to the greatest extent possible. The stack contains five records: EDC, Dome Fuji (DF), Talos Dome (TD), EPICA Dronning Maud Land (EDML) and WAIS Divide (WD). The drilling site locations are shown in Fig. 1. Volcanic ties between WD and EDC, WD and TD, and WD and EDML are developed in Buizert et al. (2018); previously published volcanic ties were used between EDC and DF (Fujita et al.2015), placing all of the records on the WD2014 chronology (Buizert et al.2015). Notably, the Vostok record, included in the stack used by Parrenin et al. (2013), is excluded from the Buizert et al. (2018) stack: it contains additional chronological uncertainty as it is derived using records from two drilling sites. We take the quadratic sum of the synchronization error from Buizert et al. (2018) and the Δage uncertainty from WD2014 to calculate the relative chronological error between ATS3 and the WD CO2 record (Fig. 2).

Figure 2Relative chronological uncertainty between ATS3 and the WD CO2 record (red line) calculated as the quadratic sum of the synchronization error from Buizert et al. (2018) (black line) and the Δage uncertainty from WD2014 (dashed blue line).


2.2CO2 and air chronology

We use atmospheric CO2 data from the WD ice core (Marcott et al.2014) which consist of 1030 measurements at 320 depths that correspond to ages between 23 000 and 9000 years BP with a median resolution of 25 years. At WD, the age offset between the ice and air (trapped much later) at a given depth, Δage, is calculated using a firn densification model, which is constrained using nitrogen-15 data, a proxy for firn column thickness (Buizert et al.2015). Δage ranges from 500±100 years at the last glacial maximum, to 200±30 years during the Holocene. Δage uncertainty is added to cumulative layer counting uncertainty to determine the total uncertainty in the air chronology.

2.3 Identifying changes in trends

We identify likely change points using piecewise linear functions. Residuals are calculated between the raw data and linear functions with a fixed number of stochastically proposed change points, which are free to explore the entire temporal range of the time series (similarly to Parrenin et al.2013). These residuals are summed to form a cost function, which allows us to perform a Bayesian analysis of the probable timing of change points. At the base of our method is a parallelized Metropolis–Hastings (MH) procedure (Goodman and Weare2010; Foreman-Mackey et al.2013). Therefore, we do not present a single “best fit” but rather analyze the ensemble of fits accepted by the routine. We plot two histograms: an upward-oriented histogram for concave-up change points and a downward-oriented histogram for concave-down change points. We use these histograms as probabilistic locators of changes in slope (Fig. 3).

The change point representations of the ATS3 and CO2 time series are composed of a set of n specified change points Xi=(xi,yi)|i=1,,n. We denote the vector of m time series observations o at time t Ol=(tl,ol)|l=1,,m, and the scalar residual term J between observations and the linear interpolation between change points fy:

(1) J ( X i ) = R T C - 1 R ; r l = f y ( t l ) - o l l ,

where R is the vector of residuals at each data point with components rl, and C is the covariance matrix of the residuals. The ATS3 series contains 1412 data points and the WD CO2 series contains 320, each of which is considered in the residuals.

We fix x0=t0 and xn=tl, i.e., the x values of the first and last change points are fixed to the first and last x values of the observation vector, with the y values allowed to vary. The remaining points are allowed to vary freely in both dimensions.

2.3.1 Estimating the covariance matrix C: treating uncertainty and noise

Our method fits time series with piecewise linear functions, and the residual vector thus accounts for any variability that cannot be represented by these fits. Paleoclimate time series, like the CO2 and ATS3 series used here, typically contain autocorrelated noise (see Mudelsee (2002), for example) which cannot be accurately represented by a piecewise linear function. Weighting the residuals of a cost-function-based formulation by a properly estimated inverse covariance matrix ensures that this autocorrelated noise is not overfitted, and can improve the balance of precision and accuracy of the fits.

Our time series contain two potential sources of uncertainty: measurement or observational uncertainty, related with the creation of the data series, and modeling uncertainty, related to the formulation of the fitting function. We formulate a separate covariance matrix to account for each source of uncertainty. These matrices are then summed to form C. We assume the measurement uncertainty to be uncorrelated in time (i.e., a white noise process). Thus, the associated covariance matrix Cmeas is diagonal and the diagonal elements Cjj are each equal to the variance of observation oj and σj2, as estimated during the measurement process.

The covariance matrix of the modeling uncertainty, which we denote as Cmod, is more complicated since the residual vector contains any autocorrelated noise in the time series that is not accounted for by the piecewise linear fits. Additionally, the time series contain outliers with respect to these linear fits and these can impact any nonrobust estimate of covariance. Finally, an initial idea of the model must be used to calculate residuals, and thus estimate their covariance. These challenges can be circumvented when data resolution is low enough to assume that residuals are uncorrelated, as in Parrenin et al. (2013); however, including the covariance matrix allows us to make use of noisy high-resolution data.

We arrive at an initial model by running a MH simulation in which C is assumed equal to the identity matrix, and select the best fit of this run. Note that Cmeas is not taken into account at this point since we require an independent estimate of Cmod. At this point, covariance could be estimated directly but tests indicated that this method was not robust, making the covariance matrix estimate sensitive to outliers and to the initial model fit. Our CO2 data are unevenly spaced in time and developing a covariance matrix using the traditional covariance estimator would require some form of interpolation, which can introduce substantial error.

The residuals with respect to the initial model are instead used to fit an AR(1) model (Robinson1977; Mudelsee2002) which treats the autocorrelation between a pair of residuals ri and ri−1 as a function of the separation between the two data points in time, ti-ti-1. The Robinson (1977) and Mudelsee (2002) model is expressed as follows:

(2) r i = r i - 1 a t i - t i - 1 ,

where the constant a determines the correlation between two residuals separated by ti-ti-1 units of time, and minimizing the loss function:

(3) S ( a ) = i = 1 n { r i - r i - 1 a t i - t i - 1 }

allows us to estimate a. We do so using a nonlinear least-squares estimate with L1-norm regularization to provide a robust estimate (Chang and Politis2016). We test the validity of the AR(1) hypothesis by comparing ri with ri-1ati-ti-1 (the Supplement). Given that the AR(1) hypothesis cannot be rejected, we can use a to calculate the theoretical correlation between two residuals, and construct a correlation matrix K and the model covariance matrix Cmod as follows:

(4) C mod = σ mod 2 K ; K i j = a t j - t i ,

where σmod2 is the variance in the modeling error, assumed constant and estimated using a robust estimator based on the interquartile range (IQR), calculated as (IQR(R)∕1.349)2 (Ghosh2018; Silverman1986). Finally, the covariance matrix of the residuals C is calculated as

(5) C = C mod + C meas .

Rather than inverting the covariance matrix, we use Cholesky and lower–upper (LU) decompositions to solve for the cost function value J, as in Parrenin et al. (2015).

2.4 Estimating the posterior probability density

In general, the probability density of the change points cannot be assumed to follow any particular distribution, as short-timescale variations of the time series may lead to multiple modes or heavy tails, for example. Thus, stochastic methods, which are best adapted to exploring general probability distributions (for example, Tarantola2005), are suited to our problem.

To tackle the large computation time required for traditional MH sampling, we apply the ensemble sampler developed by Goodman and Weare (2010) (GW) as implemented in the Python emcee library (Foreman-Mackey et al.2013). This sampler adapts the MH algorithm so that multiple model walkers can explore the probability distribution at once, making the algorithm parallelizable. It has the advantage of being affine invariant: that is, steps are adapted to the scale of the posterior distribution in a given direction.

The final task in our piecewise linear analysis is to identify the number of change points to best represent the two series we wish to analyze. The choice should reflect our goal of accurately investigating millennial-scale variability. Further, we aim for parsimony in the representation. To best balance these two goals, we apply the Bayesian information criterion (BIC; Schwarz1978) to the number of points we allow to fit the two series.

We apply a joint BIC – normalizing the cost function for each series by its lowest value – and arrive at the conclusion that the two series are best compared by fitting eight points. The histograms created for fits of seven, eight and nine points of the two series are remarkably similar, and we assess that our choice of eight points does not add significant uncertainty to the timing of change points, with the exception of the change point at the ACR onset. We include histograms of fits between five and nine points in the Supplement. We also include change point timings and lead-lag estimates calculated using seven-point fits in the Supplement.

The most probable timings are identified by probability peaks, or modes, for a fit of n points; we analyze the n time periods with greatest contiguous cumulative probability. Thus, we analyze a coherent number of change points, and avoid setting artificial probability thresholds. We avoid comparing incoherent modes by separating changes by the sign of the change in slope of the fits. If the slope decreases at a change point, the change in slope is negative or concave-down. These changes are indicated by the downward-facing part of the histogram graphs. Note that while this part of the histogram appears “negative”, probabilities cannot be negative; and this simply indicates that the probability is for a concave-down change point. If the slope increases at a change point, the change in slope is positive or concave-up. The probabilities of these changes are indicated by the upward-facing or “positive” part of the histogram. When we calculate leads and lags, we only do so for either a region in which there is a probability peak for a concave-down change point in both series, or for a concave-up change point in both series, but we do not treat concave-up probability and concave-down probability together.

2.5 Phasing

We estimate ρleadATS3, the probability that ATS3 leads CO2 over a given interval as

(6) ρ lead ATS 3 = ( ρ x ATS 3 ρ x CO 2 ) ρ chron ,

where ρxATS3 is the probability of a change point at time x for ATS3; ρxCO2 is the probability of a change point at time x for CO2; is the cross-correlation operator, which is used to calculate the probability of the difference between two variables; and is the convolution operator, which is used to calculate the probability of the sum of two variables. ρchron is the probability distribution of the chronological uncertainty between the two records, which we take to be Gaussian centered on 0, with standard deviation σ = σchron (shown in Fig. 3). The intervals associated with each change point are given in Fig. 6.

3 Results and discussion

3.1 Change point timings

Figure 3(a) Atmospheric CO2 (black) and ATS3 (red) placed on a common timescale with the normalized histograms of probable change points (eight-point simulations, allowed to fit six points and two endpoints per series). Histograms are plotted downward-oriented when the rate of change decreases and upward-oriented when it increases (same colors, y axis not shown). Probabilities are normalized so that the integrated probability for a given histogram sums to 1. In four distinct time intervals, both series show concurrent probable change points. We also plot WD Acidity (army green) and WD CH4 (violet) series. CH4 tracks changes in Northern Hemisphere climate. CO2 modes correspond with rapid changes in CH4 at the ACR end, ACR onset, 16 ka rise and the rapid rise preceding the T1 end. (b) Chronological uncertainty, taken as the sum of the Δage uncertainties and the uncertainty estimate for our volcanic synchronization.


The change point histograms for the ATS3 and CO2 time series in Fig. 3 confirm that the millennial-scale changes in the two series were largely coherent. We focus on four major changes in trends which are common to both series: the onset of the deglaciation from 18.2 to 17.2 ka, the onset of the Antarctic Cold Reversal (ACR) at around 14.5 ka, the ACR end between 12.9 and 12.65 and the end of the deglaciation at approximately 11.5 ka. For each of these four changes, we calculate the probability of a lead or lag over the time interval that encompasses the continuous peaks in the two histograms. Two CO2 change points, one centered at approximately 16 ka, and one just before the ACR onset at 14.4 ka, do not have high-probability counterparts in the ATS3 series. A low-probability change point for the temperature series after the ACR onset, centered at 14 ka, does not have a counterpart in the CO2 series.

The deglaciation onset begins with a large, positive change point mode for Antarctic temperature, centered around 18.1 ka. The corresponding change point for the CO2 series is centered around 17.6 ka.

The CO2 rise peaks at around 16 ka, identified by a downward-oriented probability peak, which has no counterpart in the temperature histogram. This peak is followed by a brief plateau in CO2 concentrations, before a gradual, accelerating resumption of the increase.

At the onset of the Antarctic Cold Reversal, an upward-facing CO2 change point at 14.4 ka is followed by a broad, downward-facing CO2 change point which peaks at around 14.25 ka. The first peak appears to reflect a centennial-scale change (identified by Marcott et al.2014) and the broadness of the second peak reflects further methodological uncertainty with respect to the timing of the millennial-scale change in CO2. An unambiguous negative temperature change also occurs at around 14.3 ka, roughly concurrent with the downward CO2 change point. Antarctic temperature began to descend after the ACR onset, and it is worth mentioning the low-probability concave-up change point mode centered on 13.9 ka, particularly because this point is much more probable for some of the individual isotopic records. No corresponding change point is detected for CO2.

The ACR termination is represented by large probability peaks in both series. An increase in CO2 began at the peak occurring around 12.85 ka, while the ATS3 increase is centered at approximately 12.78 ka reaching its maximum around 12.7 ka.

The end of the deglacial warming in Antarctica is well-defined in the ATS3 series, with a large mode reaching its maximum at 11.7 ka. The corresponding CO2 mode reaches its maximum at around 11.15 ka. Visually, we might question why the CO2 change point is deemed to most probably occur at 11.15 ka, when a kink in the series appears to occur closer to 11.5 ka. It is important to note that two minor spikes in probability appear to fit the rapid rise that occurs close to 11.5 ka, and that the downward spike at the end of this rise indeed indicates that a small number of accepted iterations do indeed fit a change point here.

Figure 4(a) Savitzky–Golay filtered atmospheric CO2 (black) and ATS3 (red) placed on a common timescale with the normalized histograms of probable change points (eight-point simulations, allowed to fit six points per series). Histograms are plotted downward-oriented when the rate of change decreases and upward-oriented when it increases (same colors, y axis not shown, probabilities range from 0 (center) to 0.0024 (top and bottom)). (b) Chronological uncertainty, taken as the sum of the Δage uncertainties and the uncertainty estimate for our volcanic synchronization.


Since the resolution of the WAIS Divide dataset decreases considerably after the Holocene onset, and only three points account for the majority of the rapid rise that occurs before the Holocene onset, most of the weight is given to the obvious line beginning at the ACR end. Adding an additional point should allow us to fit this slightly better; the CO2 series is slightly better fit with nine points according to the individual BIC values, and there is more probability around the rapid rise in the nine-point fit, though the peak at 11.15 ka is still dominant. In any case, some methodological uncertainty exists regarding the location of this point, and the probability estimate is possibly biased by the quick change in resolution. Better resolution around this point will help identify the true location of the change.

As a second test of the timings of millennial-scale events, we use our method to fit filtered versions of the ATS3 and WAIS Divide CO2 data. A Savitzky–Golay filter, designed to have an approximate cutoff periodicity of 500 years, is applied to the two records. Fitting change points to these two series allows us to verify that our leads and lags are not overly influenced by submillennial scale noise in the original records.

Figure 4 shows the Savitzky–Golay filtered CO2 and ATS2 time series, and the corresponding change point histograms. The four major changes identified in both series, at the T1 onset, the ACR onset, the ACR end and T1 end, are similar in shape and center to the change points identified for the raw data. However, there are two notable differences between the two fits. First, the histograms are smoother and have broader peaks. This is not surprising, given that the Savitzky–Golay filters are designed to remove all variability with periodicities less than 500 years, whereas the covariance matrix applied to the fits of the raw data only treats an approximation of AR(1)-correlated noise. Second, the pre-ACR change in CO2 is removed from the filtered series, which is again reasonable as it appears to mark a centennial-scale event. Savitzky–Golay filtering has its own drawbacks – data reinterpolation is required, for example, and propagating measurement uncertainty becomes difficult. However, the similarity of the two results supports our fits of the raw data.

3.2 Change point timings for individual temperature records

Figure 5Atmospheric CO2 (black) and individual δ18O records placed on a common timescale, with the normalized histograms of probable change points (eight points) for each ice core used in the ATS3 stack; the locations of the drill sites are shown in the top right corner. Details of the histogram plots are as in Fig. 3. Maximum probability estimates and 68 % and 95 % probability intervals for timings of the individual records are provided in the Supplement. The δ18O records are given in per mill anomalies with respect to the last 200 years, as is the ATS3 stack.


Histograms calculated for each of the regional δ18O records are shown in Fig. 5. These histograms should still be interpreted cautiously, as additional information included in the isotopic records here assumed to represent temperature – the signal of ice sheet elevation change, for example – are not corrected for. The comparison of these histograms provides an initial, exploratory picture of potential regional differences in climate change during the last termination.

Of the four changes identified as coherent between the temperature stack and CO2, those at the deglaciation onset, the ACR end and the T1 end are expressed as probability peaks in all five records. Some ambiguity appears to exist about the timing of the ACR onset in the EDML record. It is expressed by a rather broad, low-probability mode extending between 16 and 14 ka, though a large spike at 14 ka marks the downturn seen in the other records. The ACR onset is well-defined in all of the five other records. Three of the records – TD, EDC and WD – show a marked stabilization in temperature after the ACR onset, near 13.8 ka, which appears as a region of much lower probability in the ATS3 stack.

The WAIS Divide record is, notably, the only isotopic record in the stack from the West Antarctic Ice Sheet. We could thus reasonably expect it to show considerably different trends from the other records. Indeed, changes in the WD temperature record occur at 22 and 20 ka. This earlier change in the isotopic record was identified and confirmed to indeed be a temperature signal by Cuffey et al. (2016) using a borehole temperature record, though their study places the change at 21 ka. We confirm that the onset of the deglacial temperature rise in West Antarctica likely began as much as 4 kyr before the onset of temperature rise in East Antarctica. Interestingly, the WD record also shows a temperature change point around 17.8 ka, expressed slightly later than in the other records and more synchronous with CO2. This apparent acceleration of the temperature rise is followed by a downward-facing change point not seen in any of the other records. A difference appears to exist in timing at the T1 end as well, with temperature change at WD appearing to precede the East Antarctic records and the DF temperature change, centered at 11.2 ka, occurring more synchronously with CO2.

3.3 Leads and lags

The probability densities of leads and lags at the coherent change points between ATS3 and CO2 are shown in Fig. 6. We report the central 68 % and 95 % probability intervals for each histogram. These values are grouped in Table 1.

Figure 6Probability density ρ (y axis, normalized) of an ATS lead (x-axes, in years) at each of the selected change point intervals (noted on subfigures). Negative x-axis values indicate a CO2 lead. The maximum probability lead/lag and 68 %95 % central probability intervals are indicated by the orange dot and lines on each histogram.


Parrenin et al. (2013)

Table 1Maximum probabilities and central probability intervals for leads and lags at each of the selected change point intervals. Negative x-axis values indicate a CO2 lead.

a Note the small, but distinct, probability that the CO2 change point occurs closer to 11.5 ka would indicate a lag of 174 years, with a 68 % central probability range of 65 to 280 years. b The phasing at the ACR onset is sensitive to whether seven or eight points are used. Timings with seven- and eight-point fits, and for the individual isotopic records, are made available in the Supplement.

Download Print Version | Download XLSX

ATS3 led CO2 by 570 years, (within a 68 % interval of 128 to 751 years) at the T1 onset. Given the large range of uncertainty, though, we cannot exclude the possibility of synchrony at the 95 % level, which, interestingly, appears to be the case for the Dome Fuji record. At the ACR onset, we are not able to identify a clear lead or lag. At this point, phasing is sensitive to the number of points used to make the calculation: with 7 points, we calculate a 240 year lead of ATS3, and with 8 points, we calculate a 50 year lead. In neither of these cases can we exclude synchrony within 95 % probability, and with 8 points, it is well within 68 %. At the ACR end, CO2 led ATS3, by 65 ± years within a 68 % range of 211 years to −117 years (a temperature lead) and so again, the possibility of synchrony cannot be excluded within 68 % probability.

At the T1 end, a CO2 lag is certain. Calculating the phasing between 12.0 and 11.0 ka, we obtain an ATS3 lead of 532 years, with a 68 % probability range of 337 to 629 years. This estimate is complicated, though, if we consider the small possibility that the true CO2 change point occurs closer to 11.5 ka, at the end of the rapid rise. In this case, the phasing is reduced to 174 years (68 % central probability range of 65 to 280 years) and synchrony is within the 95 % central probability interval (−71 to 411 years).

3.4 Discussion

Our results refine and complicate the timings, leads and lags identified by the most recent comparable studies (Pedro et al.2012; Parrenin et al.2013). We identify a CO2 change point not treated in these studies at 16 ka and one before the ACR onset, associated with the centennial-scale rapid rises identified by Marcott et al. (2014). We also treat regional isotopic records and identify a change point occurring at 13.9 ka in three of the records. None of these change points have a marked counterpart in the other series.

During the major multi-millennial-scale changes which occur at T1 onset and T1 end, Antarctic temperature likely led CO2 by several centuries. However, during the complex centennial-scale change at the ACR onset, we cannot calculate a clear lead of either ATS or CO2; and at the end of the ACR, CO2 leads temperature. Further, we neither identify a temperature analog for the CO2 change at 16 ka nor an analog in CO2 of the temperature stabilization in ATS3 after the ACR onset, itself not present in all of the regional δ18O records, indicating at least some degree of decoupling during these changes. Additionally, the CO2 changes at the ACR onset and T1 end are overlaid with centennial-scale substructures. Finally, synchrony is within the 2σ uncertainty range for each phasing, with the exception of the T1 end.

The changes in CO2 occurring at the ACR onset, ACR end and the T1 end have been identified to correspond with changes in CH4 (Marcott et al.2014), which are thought to originate in tropical wetland sources (Chappellaz et al.1997; Fischer et al.2008; Petrenko et al.2009) and are indicative of Northern Hemisphere and low-latitude temperature changes during the deglaciation (Shakun et al.2012). Indeed, the CO2 modes appear to demarcate the rapid changes in the WD CH4 record, shown in Fig. 3.

The beginning of a gradual rise in CH4 at around 18 ka appears to be near-synchronous with the T1 onset rise in Antarctic temperature. This rise is not seen in Greenland paleotemperature records, where it may have been masked by AMOC-driven wintertime cooling (Buizert et al.2017) but it also appears in proxy temperature stacks spanning both the Northern Hemisphere and Southern Hemisphere 0 to 30 latitude bands (Shakun et al.2012).

Tephras from Mt. Takahe, a stratovolcano located in West Antarctica, have been detected in Antarctic ice cores during a 192-year interval around 17.7 ka. It has been postulated that this eruption may have provoked changes to large-scale SH circulation via ozone depletion, possibly triggering the transition between the gradual SH temperature rise beginning well before 18 ka and the more rapid rise marking the deglaciation (McConnell et al.2017). The CO2 mode we find at the deglaciation is coeval with this event within the range of dating uncertainty (Fig. 3), and CH4 visually appears to accelerate concurrently. However, the cumulative probability of the ATS3 change point is much greater before 17.7 ka than after (approximately 80 % of the probability density occurs before, see Fig. 3); hence, our results do not support the proposed volcanic forcing by McConnell et al. (2017) of the temperature change.

Though the T1 onset and the ACR end are both thought to originate in AMOC reductions (Marcott et al.2014), our results allow for the CO2–ATS3 phasing to be different during the two events, with the maximum probabilities reversed in directionality (i.e. with temperature leading at T1 and CO2 leading at the ACR end, though zero phasing is within 95 % error). Though CH4 appears to change alongside CO2 during both intervals, the phasing between CO2 and ATS3 are opposite in direction and different in slope. This hints at a complex coupling, depending on conditions defined by multiple variables and mechanisms between CO2 and Antarctic temperature. Bauska et al. (2016), for example, hypothesize that an earlier rise of CO2 at 12.9 ka, driven by land carbon loss or SH westerly winds, might have been superimposed on the millennial-scale trend.

The apparent decoupling between CO2 and ATS3 at 16 ka also merits further discussion. We do not detect a change point for any of the isotopic records at 16 ka but the EDML record contains extremely broad uncertainty associated with the ACR onset peak, stretching to 16 ka, which indicates that this portion of the EDML time series is indeed notably different in shape from the other records, even if a clear signal is not identified at 16 ka by our method. EDML indeed appears to record changes in AMOC differently than the other isotopic records (Buizert et al.2018; Landais et al.2018). CO2 and ATS3 are similarly apparently decoupled at the temperature change point centered at 14 ka and this point could be indicative of variability specific to the Pacific and eastern Indian Ocean sectors, as it is present only in the TALOS Dome and EDC records, and slightly later, around 13.7 ka, in the WAIS Divide record, indicating a cooling trend after the ACR onset which is not clear in the DF or EDML series.

Within the range of uncertainty, the mean of our lead–lag estimates is consistent with the boundaries proposed by Pedro et al. (2012). Our results are consistent with those of Parrenin et al. (2013) for three out of the four change points addressed, but differ considerably at the T1 onset.

The considerable difference at T1 between our result and that of Parrenin et al. (2013) is most likely due to the much higher resolution of the WD CO2 time series. It is also possible that the result of Parrenin et al. (2013) was limited to a local probability maximum of this change point in the CO2 series. The addition of the WD paleotemperature record and removal of the Vostok record from ATS3, the updated atmospheric CO2 dataset, and our more generalized methodology are all, in part, responsible for the differences in computed time delays (the Supplement). This testifies to the importance of data resolution, methodological development, and chronological accuracy in the determination of leads and lags.

4 Conclusions

Our study is a follow-up of the studies by Pedro et al. (2012) and Parrenin et al. (2013) on the leads and lags between atmospheric CO2 and Antarctic temperature during the last deglacial warming. We refine the results of these studies by using the high-resolution CO2 record from WD; using Δage computed on WD; using a new Antarctic temperature stack composed of 5 volcanically synchronized ice core isotope records, developed by Buizert et al. (2018); and using a more precise and complete probabilistic estimate to determine change points. Our methodology detects four major common break points in both time series. The phasing between CO2 and Antarctic climate is small but variable, with phasing ranging from a centennial-scale CO2 lead, to synchrony, to a multicentennial-scale lead of Antarctic climate. This variability in phasing indicates that the mechanisms of coupling are complex. We propose three possibilities: (i) the mechanisms by which CO2 and Antarctic temperature were coupled were consistent through the deglaciation, but can be modulated by external forcings or background conditions that impact heat transfer and oceanic circulation (and hence CO2 release); (ii) these mechanisms can be modulated by internal feedbacks that change the response timings of the two series; and/or (iii) multiple, distinct mechanisms might have provoked similar responses in both series, but with accordingly different lags.

We also explore the hypothesis of regional differences in temperature change in Antarctica. Though the use of individual isotopic temperature records is complicated by influences other than regional temperature, including localized variations in source temperature and ice sheet elevation change we confirm that the deglacial temperature rise did not occur homogeneously across the Antarctic continent, with significant differences existing between the WAIS Divide and East Antarctic records at the onset of the termination and smaller potential differences occurring between the East Antarctic records, including a considerably later end of the deglacial warming in the Dome Fuji record.

Hypotheses of relationships between these events should now be reinvestigated with modeling studies. The relationship between CO2 and Antarctic temperature on longer timescales and during other periods of rapid climate change is also of interest. Additional high-resolution West Antarctic paleotemperature records would allow for a robust investigation of regional differences between West Antarctica and East Antarctica and our analysis at the T1 end could be improved with continued high-resolution CO2 measurements through the beginning of the Holocene. Finally, the continued measurement of high-resolution ice core CO2 records is essential to understand the relationship between CO2 and global and regional temperature during the last 800 000 years.

Code and data availability

The code and data series used in this article are available at (Chowdhry Beeman, 2019). The CO2, ATS3 and individual isotopic series, as well as the chronological uncertainty, are included as .txt files. The original CO2 data from Marcott et al. (2014) are available at The synchronization tie points and original ATS3 stack are available in the Supplementary Data spreadsheet of Buizert et al. (2018) at


The supplement related to this article is available online at:

Author contributions

JCB and LG contributed equally to this manuscript, and should both be considered lead authors. LG, FP, DR and JCB designed the research project. JCB and FP designed the change point identification method. JCB and LG led the writing of the manuscript; and all authors contributed to analysis and discussion of the phasing and timings, with TJF, CB and EB in particular providing expertise on the ATS3 stack, the WD CO2 record and their climatic context; and DR and FP on the lead–lag problem and its history.

Competing interests

The authors declare that they have no conflict of interest.


We thank Michael Sigl, Jinhwa Shin, Emmanuel Witrant and Amaelle Landais for their support and great help discussing this work and Mirko Severi for his EDC data and support with the volcanic synchronization. This work is supported by the Fondation Ars et Cuttoli and by the LEFE IceChrono and CO2Role projects.

Financial support

This research has been supported by the Centre National de la Recherche Scientifique (grant LEFE IceChrono and LEFE CO2Role) and the Fondation Ars et Cuttoli (grant CO2Role).

Review statement

This paper was edited by Hubertus Fischer and reviewed by Joel Pedro and one anonymous referee.


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. a, b

Anklin, M., Barnola, J.-M., Schwander, J., Stauffer, B., and Raynaud, D.: Processes affecting the CO2 concentrations measured in Greenland ice, Tellus B, 47, 461–470, 1995. a

Barnola, J.-M., Pimienta, P., Raynaud, D., and Korotkevich, Y. S.: CO2-climate relationship as deduced from the Vostok ice core: a re-examination based on new measurements and on a re-evaluation of the air dating, Tellus B, 43, 83–90, 1991. a

Bauska, T. K., Baggenstos, D., Brook, E. J., Mix, A. C., Marcott, S. A., Petrenko, V. V., Schaefer, H., Severinghaus, J. P., and Lee, J. E.: Carbon isotopes characterize rapid changes in atmospheric carbon dioxide during the last deglaciation, P. Natl. Acad. Sci. USA, 113, 3465–3470, 2016. a, b

Berger, A.: Long-term variations of daily insolation and Quaternary climatic changes, J. Atmos. Sci., 35, 2362–2367, 1978. a

Buizert, C. and Severinghaus, J. P.: Dispersion in deep polar firn driven by synoptic-scale surface pressure variability, The Cryosphere, 10, 2099–2111,, 2016. a, b

Buizert, C., Cuffey, K. M., Severinghaus, J. P., Baggenstos, D., Fudge, T. J., Steig, E. J., Markle, B. R., Winstrup, M., Rhodes, R. H., Brook, E. J., Sowers, T. A., Clow, G. D., Cheng, H., Edwards, R. L., Sigl, M., McConnell, J. R., and Taylor, K. C.: The WAIS Divide deep ice core WD2014 chronology – Part 1: Methane synchronization (68–31 ka BP) and the gas age–ice age difference, Clim. Past, 11, 153–173,, 2015. a, b, c

Buizert, C., Keisling, B. A., Box, J. E., He, F., Carlson, A. E., Sinclair, G., and DeConto, R. M.: Greenland-Wide Seasonal Temperatures During the Last Deglaciation, Geophys. Res. Lett., 45, 1–10,, 2017. a

Buizert, C., Sigl, M., Severi, M., Markle, B. R., Wettstein, J. J., McConnell, J. R., Pedro, J. B., Sodemann, H., Goto-Azuma, K., Kawamura, K., et al.: Abrupt ice-age shifts in southern westerly winds and Antarctic climate forced from the north, Nature, 563, 681–685, 2018. a, b, c, d, e, f, g, h

Burke, A. and Robinson, L. F.: The Southern Ocean's role in carbon exchange during the last deglaciation, Science, 335, 557–561, 2012. a

Caillon, N., Severinghaus, J. P., Jouzel, J., Barnola, J.-M., Kang, J., and Lipenkov, V. Y.: Timing of atmospheric CO2 and Antarctic temperature changes across Termination III, Science, 299, 1728–1731, 2003. a

Chang, C. C. and Politis, D. N.: Robust autocorrelation estimation, J. Comput. Graph. Stat., 25, 144–166, 2016. a

Chappellaz, J., Blunier, T., Kints, S., Dällenbach, A., Barnola, J.-M., Schwander, J., Raynaud, D., and Stauffer, B.: Changes in the atmospheric CH4 gradient between Greenland and Antarctica during the Holocene, J. Geophys. Res.-Atmos., 102, 15 987–15 997, 1997. a

Chowdhry Beeman, J.: Jai-Chowdhry/LinearFit-2.0-beta 2.0-CP (Version 2.0-CP), Zenodo, available at:, last access: 7 May 2019. 

Cuffey, K. M., Clow, G. D., Steig, E. J., Buizert, C., Fudge, T., Koutnik, M., Waddington, E. D., Alley, R. B., and Severinghaus, J. P.: Deglacial temperature history of West Antarctica, P. Natl. Acad. Sci. USA, 113, 14249–14254, 2016. a, b, c

Fischer, H., Behrens, M., Bock, M., Richter, U., Schmitt, J., Loulergue, L., Chappellaz, J., Spahni, R., Blunier, T., Leuenberger, M., and Stocker, T. F.: Changing boreal methane sources and constant biomass burning during the last termination, Nature, 452, 864–867, 2008. a

Foreman-Mackey, D., Hogg, D. W., Lang, D., and Goodman, J.: Emcee: The MCMC Hammer, Publ. Astron. Soc. Pac., 125, 306–312,, 2013. a, b

Fujita, S., Parrenin, F., Severi, M., Motoyama, H., and Wolff, E. W.: Volcanic synchronization of Dome Fuji and Dome C Antarctic deep ice cores over the past 216 kyr, Clim. Past, 11, 1395–1416,, 2015. a

Ghosh, S.: Kernel smoothing: Principles, methods and applications, John Wiley & Sons, Hoboken, NJ, USA, 2018. a

Goodman, J. and Weare, J.: Ensemble samplers with affine invariance, Comm. App. Math. Com. Sc., 5, 65–80, 2010. a, b

Hays, J. D., Imbrie, J., and Shackleton, N. J.: Variations in the Earth's orbit: pacemaker of the ice ages, Science, 194, 1121–1132, 1976. a

Johnsen, S., Dansgaard, W., and White, J.: The origin of Arctic precipitation under present and glacial conditions, Tellus B, 41, 452–468, 1989. a

Jones, J. M., Gille, S. T., Goosse, H., Abram, N. J., Canziani, P. O., Charman, D. J., Clem, K. R., Crosta, X., De Lavergne, C., Eisenman, I., England, M. H., Fogt, R. L., Frankcombe, L. M., Marshall, G. J., Masson-Delmotte, V., Morrison, A. K., Orsi, A. J., Raphael, M. N., Renwick, J. A., Schneider, D. P., Simpkins, G. R., Steig, E. J., Stenni, B., Swingedouw, D., and Vance, T. R.: Assessing recent trends in high-latitude Southern Hemisphere surface climate, Nat. Clim. Change, 6, 91–926, 2016. a

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, J., 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. a, b, c

Kawamura, K., Parrenin, F., Lisiecki, L., Uemura, R., Vimeux, F., Severinghaus, J. P., Hutterli, M. A., Nakazawa, T., Aoki, S., Jouzel, J., Raymo, M. E., Matsumoto, K., Nakata, H., Motoyama, H., Fujita, S., Goto-Azuma, K., Fujii, Y., and Watanabe, O.: Northern Hemisphere forcing of climatic cycles in Antarctica over the past 360,000 years, Nature, 448, 912–916, 2007. a

Landais, A., Capron, E., Masson-Delmotte, V., Toucanne, S., Rhodes, R., Popp, T., Vinther, B., Minster, B., and Prié, F.: Ice core evidence for decoupling between midlatitude atmospheric water cycle and Greenland temperature during the last deglaciation, Clim. Past, 14, 1405–1415,, 2018. a

Lisiecki, L. and Raymo, M.: A Plio-Pleistocene stack of 57 globally distributed benthic 0126-026180 records, Paleoceanography, 20, 1–17, 2005. a

Lorius, C. and Merlivat, L.: Distribution of mean surface stable isotopes values in East Antarctica; observed changes with depth in coastal area, Tech. rep., CEA Centre d'Etudes Nucleaires de Saclay, Gif-Sur-Yvette, France, 1975. a

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. a

Lüthi, 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. a

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

Masson-Delmotte, V., Schulz, M., Abe-Ouchi, A., Beer, J., Ganopolski, A., Gonzaález Rouco, J., Jansen, E., Lambeck, K., Luterbacher, J., Naish, T., Osborn, T., Otto-Bliesner, B., Quinn, T., Ramesh, R., Rojas, M., Shao, X., and Timmermann, A.: Information from Paleoclimate Archives, in: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T., Qin, D., Plattner, G.-K., Tignor, M., Allen, S., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P., book section 5, 383–464, Cambridge University Press, Cambridge, UK and New York, NY, USA,, 2013. a

McConnell, J. R., Burke, A., Dunbar, N. W., Köhler, P., Thomas, J. L., Arienzo, M. M., Chellman, N. J., Maselli, O. J., Sigl, M., Adkins, J. F., Baggenstos, D., Burkhart, J. F., Brook, E. J., Buizert, C., Cole-Dai, J., Fudge, T. J., Knorr, G.,Graf, H.-F., Grieman, M. M., Iverson, N., McGwire, K. C., Mulvaney, R., Paris, G., Rhodes, R. H., Saltzman, E. S., Severinghaus, J. P., Steffensen, J. P., Taylor, K. C., and Winckler, G.: Synchronous volcanic eruptions and abrupt climate change ∼17.7 ka plausibly linked by stratospheric ozone depletion, P. Natl. Acad. Sci. USA, 114, 10035–10040, 2017. a, b

Monnin, E., Indermühle, A., Dällenbach, A., Flückiger, J., Stauffer, B., Stocker, T. F., Raynaud, D., and Barnola, J.-M.: Atmospheric CO2 concentrations over the last glacial termination, Science, 291, 112–114, 2001. a

Mudelsee, M.: TAUEST: A computer program for estimating persistence in unevenly spaced weather/climate time series, Comput. Geosci., 28, 69–72, 2002. a, b, c

NorthGRIP Project Members: High-resolution record of Northern Hemisphere climate extending into the last interglacial period, Nature, 431, 147–151, 2004. a

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. a

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. a, b, c, d, e, f, g, h, i, j, k, l

Parrenin, F., Bazin, L., Capron, E., Landais, A., Lemieux-Dudon, B., and Masson-Delmotte, V.: IceChrono1: a probabilistic model to compute a common and optimal chronology for several ice cores, Geosci. Model Dev., 8, 1473–1492,, 2015. a

Pedro, J. B., Rasmussen, S. O., and van Ommen, T. D.: Tightened constraints on the time-lag between Antarctic temperature and CO2 during the last deglaciation, Clim. Past, 8, 1213–1221,, 2012. a, b, c, d, e, f

Pedro, J. B., Bostock, H. C., Bitz, C. M., He, F., Vandergoes, M. J., Steig, E. J., Chase, B. M., Krause, C. E., Rasmussen, S. O., Markle, B. R., and Cortese, G.: The spatial extent and dynamics of the Antarctic Cold Reversal, Nat. Geosci., 9, 51, 2016. a

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. a

Petrenko, V. V., Smith, A. M., Brook, E. J., Lowe, D., Riedel, K., Brailsford, G., Hua, Q., Schaefer, H., Reeh, N., Weiss, R. F., Etheridge, D., and Severinghaus, J. P.: 14CH4 measurements in Greenland ice: investigating last glacial termination CH4 sources, Science, 324, 506–508, 2009. a

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

Raynaud, D. and Siegenthaler, U.: Role of trace gases: the problem of lead and lag, Global changes in the perspective of the past, edited by: Eddy, J. A., and Oeschger, H., Wiley , Chichester, UK, 173–188, 1993. a

Raynaud, D., Jouzel, J., Barnola, J., Chappellaz, J., Delmas, R., and Lorius, C.: The ice record of greenhouse gases, Science, 259, 926–926, 1993. a

Robinson, P.: Estimation of a time series model from unequally spaced data, Stoch. Proc. Appl., 6, 9–24, 1977. a, b

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

Schwarz, G.: Estimating the dimension of a model, Ann. Stat., 6, 461–464, 1978. a

Shakun, J. D., Clark, P. U., He, F., Marcott, S. A., Mix, A. C., Liu, Z., Otto-Bliesner, B., Schmittner, A., and Bard, E.: Global warming preceded by increasing carbon dioxide concentrations during the last deglaciation, Nature, 484, 49–54, 2012.  a, b, c, d

Silverman, B. W.: Density estimation for statistics and data analysis, Chapman and Hall, London, UK and Ney York, USA, 1986. a

Skinner, L., Fallon, S., Waelbroeck, C., Michel, E., and Barker, S.: Ventilation of the deep Southern Ocean and deglacial CO2 rise, Science, 328, 1147–1151, 2010. a

Sowers, T., Bender, M., Raynaud, D., and Korotkevich, Y. S.: δ15N of N2 in air trapped in polar ice: A tracer of gas transport in the firn and a possible constraint on ice age-gas age differences, J. Geophys. Res.-Atmos., 97, 15683–1697, 1992. a

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

Tarantola, A.: Inverse problem theory, SIAM, Philadelphia, USA, 2005. a

WAIS Divide Project Members: Onset of deglacial warming in West Antarctica driven by local orbital forcing, Nature, 500, 440–444, 2013. a, b

Williams, D., Peck, J., Karabanov, E., Prokopenko, A., Kravchinsky, V., King, J., and Kuzmin, M.: Lake Baikal record of continental climate response to orbital insolation during the past 5 million years, Science, 278, 1114–1117, 1997. a

Short summary
Atmospheric CO2 was likely an important amplifier of global-scale orbitally-driven warming during the last deglaciation. However, the mechanisms responsible for the rise in CO2, and the coherent rise in Antarctic isotopic temperature records, are under debate. Using a stochastic method, we detect variable lags between coherent changes in Antarctic temperature and CO2. This implies that the climate mechanisms linking the two records changed or experienced modulations during the deglaciation.