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

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 CO₂ as recorded in Antarctic ice cores. In addition to five regional records, we compare a δ¹⁸O stack, representing Antarctic climate variations with the high-resolution robustly dated WAIS Divide CO₂ record (West Antarctic Ice Sheet). We assess the CO₂ 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 CO₂ 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 CO₂ 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 (CO₂ 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 CO₂ 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 Raymo, 2005; 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; Berger, 1978; 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 CO₂ and continental-ice-surface–albedo changes. However, the mechanisms that control the CO₂ rise are still a matter of debate. Accordingly, reconstructing the phase relationship (leads and lags) between climate variables and CO₂ during the last termination has become important, and has a substantial history in ice core research (Barnola et al., 1991; Raynaud and Siegenthaler, 1993; Caillon et al., 2003; Pedro et al., 2012; Parrenin et al., 2013).

Global temperature has been shown to lag CO₂ on average during T1 (Shakun et al., 2012), supporting the importance of CO₂ as an amplifier of orbitally driven global-scale warming. But Antarctic temperature and CO₂ concentrations changed much more coherently as T1 progressed. Indeed, midway through the glacial–interglacial transition, Antarctic warming and CO₂ 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 (B∕A). 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 B∕A, 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 CO₂ increases (Anderson et al., 2009; Burke and Robinson, 2012; 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 Johnsen, 2003; 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 CO₂ 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 Members, 2004; Jouzel et al., 2007) due to the temperature-dependent fractionation of water isotopes (Lorius and Merlivat, 1975; 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 Severinghaus, 2016). 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 CO₂ 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 CO₂ (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 CO₂ began to increase (Monnin et al., 2001), a result that was sometimes misinterpreted to mean that CO₂ 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 CO₂ 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 CO₂ 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 δ¹⁵N of N₂ as a proxy of the DZ height, assuming that the height of the CZ was negligible during the study period. CO₂ and Antarctic temperature were found to be in phase at the beginning of TI ($-\mathrm{10}\pm \mathrm{160}$ years) and at the end of the ACR period ($-\mathrm{60}\pm \mathrm{120}$ years), but CO₂ 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 δ¹⁵N reflects DZ height is imperfect as it may underestimate the DZ height for sites with strong barometric pumping and layering (Buizert and Severinghaus, 2016).

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

A new CO₂ 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 Members, 2013; 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 Members, 2013).

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 CO₂. 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 CO₂ 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 δ¹⁸O 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 CO₂ record (Fig. 2).

Figure 2 Relative chronological uncertainty between ATS3 and the WD CO₂ 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.2 CO₂ and air chronology

We use atmospheric CO₂ 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 Weare, 2010; 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 CO₂ time series are composed of a set of n specified change points ${\mathit{X}}_i=(x_i,y_i)|i=\mathrm{1},\mathrm{\dots },n$. We denote the vector of m time series observations o at time t ${\mathit{O}}_l=(t_l,o_l)|l=\mathrm{1},\mathrm{\dots },m$, and the scalar residual term J between observations and the linear interpolation between change points f_y:

$$\begin{array}{}\text{(1)}& J\left({\mathit{X}}_i\right)={\mathit{R}}^T{\mathbf{C}}^{-\mathrm{1}}\mathit{R};{\mathit{r}}_l={\left(f_y\left(t_l\right)-o_l\right)}_l,\end{array}$$

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

We fix x₀=t₀ and x_n=t_l, 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 CO₂ 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 C_meas is diagonal and the diagonal elements C_jj are each equal to the variance of observation o_j and ${\mathit{\sigma }}_j^{\mathrm{2}}$, as estimated during the measurement process.

The covariance matrix of the modeling uncertainty, which we denote as C_mod, 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 C_meas is not taken into account at this point since we require an independent estimate of C_mod. 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 CO₂ 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 (Robinson, 1977; Mudelsee, 2002) which treats the autocorrelation between a pair of residuals r_i and r_i−1 as a function of the separation between the two data points in time, $t_i-t_{i-\mathrm{1}}$. The Robinson (1977) and Mudelsee (2002) model is expressed as follows:

$$\begin{array}{}\text{(2)}& r_i=r_{i-\mathrm{1}}\cdot a^{t_i-t_{i-\mathrm{1}}},\end{array}$$

where the constant a determines the correlation between two residuals separated by $t_i-t_{i-\mathrm{1}}$ units of time, and minimizing the loss function:

$$\begin{array}{}\text{(3)}& S\left(a\right)=\sum _{i=\mathrm{1}}^n\mathit{\{}{r}_i-r_{i-\mathrm{1}}\cdot a^{t_i-t_{i-\mathrm{1}}}\mathit{\}}\end{array}$$

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 Politis, 2016). We test the validity of the AR(1) hypothesis by comparing r_i with $r_{i-\mathrm{1}}\cdot a^{t_i-t_{i-\mathrm{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 C_mod as follows:

$$\begin{array}{}\text{(4)}& {\mathbf{C}}_{\mathrm{mod}}={\mathit{\sigma }}_{\mathrm{mod}}^{\mathrm{2}}\mathbf{K};{\mathbf{K}}_{ij}=a^{t_j-t_i},\end{array}$$

where ${\mathit{\sigma }}_{\mathrm{mod}}^{\mathrm{2}}$ 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)² (Ghosh, 2018; Silverman, 1986). Finally, the covariance matrix of the residuals C is calculated as

$$\begin{array}{}\text{(5)}& \mathbf{C}={\mathbf{C}}_{\mathrm{mod}}+{\mathbf{C}}_{\mathrm{meas}}.\end{array}$$

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, Tarantola, 2005), 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; Schwarz, 1978) 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 ${\mathit{\rho }}_{\mathrm{lead}}^{\mathrm{ATS}\mathrm{3}}$, the probability that ATS3 leads CO₂ over a given interval as

$$\begin{array}{}\text{(6)}& {\mathit{\rho }}_{\mathrm{lead}}^{\mathrm{ATS}\mathrm{3}}=({\mathit{\rho }}_x^{\mathrm{ATS}\mathrm{3}}\circ {\mathit{\rho }}_x^{{\mathrm{CO}}_{\mathrm{2}}})\star {\mathit{\rho }}^{\mathrm{chron}},\end{array}$$

where ${\mathit{\rho }}_x^{\mathrm{ATS}\mathrm{3}}$ is the probability of a change point at time x for ATS3; ${\mathit{\rho }}_x^{{\mathrm{CO}}_{\mathrm{2}}}$ is the probability of a change point at time x for CO₂; ∘ 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 CO₂ (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 CH₄ (violet) series. CH₄ tracks changes in Northern Hemisphere climate. CO₂ modes correspond with rapid changes in CH₄ 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 CO₂ 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 CO₂ 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 CO₂ 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 CO₂ series is centered around 17.6 ka.

The CO₂ 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 CO₂ concentrations, before a gradual, accelerating resumption of the increase.

At the onset of the Antarctic Cold Reversal, an upward-facing CO₂ change point at 14.4 ka is followed by a broad, downward-facing CO₂ 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 CO₂. An unambiguous negative temperature change also occurs at around 14.3 ka, roughly concurrent with the downward CO₂ 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 CO₂.

The ACR termination is represented by large probability peaks in both series. An increase in CO₂ 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 CO₂ mode reaches its maximum at around 11.15 ka. Visually, we might question why the CO₂ 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 CO₂ (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 CO₂ 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 CO₂ 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 CO₂ 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 CO₂ 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 5 Atmospheric CO₂ (black) and individual δ¹⁸O 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 δ¹⁸O 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 δ¹⁸O 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 CO₂, 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 CO₂. 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 CO₂.

3.3 Leads and lags

The probability densities of leads and lags at the coherent change points between ATS3 and CO₂ 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 6 Probability 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 CO₂ 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 1 Maximum probabilities and central probability intervals for leads and lags at each of the selected change point intervals. Negative x-axis values indicate a CO₂ lead.

^a Note the small, but distinct, probability that the CO₂ 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.

ATS3 led CO₂ 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, CO₂ 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 CO₂ 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 CO₂ 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 CO₂ 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 CO₂ by several centuries. However, during the complex centennial-scale change at the ACR onset, we cannot calculate a clear lead of either ATS or CO₂; and at the end of the ACR, CO₂ leads temperature. Further, we neither identify a temperature analog for the CO₂ change at 16 ka nor an analog in CO₂ of the temperature stabilization in ATS3 after the ACR onset, itself not present in all of the regional δ¹⁸O records, indicating at least some degree of decoupling during these changes. Additionally, the CO₂ 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 CO₂ occurring at the ACR onset, ACR end and the T1 end have been identified to correspond with changes in CH₄ (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 CO₂ modes appear to demarcate the rapid changes in the WD CH₄ record, shown in Fig. 3.

The beginning of a gradual rise in CH₄ 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 CO₂ mode we find at the deglaciation is coeval with this event within the range of dating uncertainty (Fig. 3), and CH₄ 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 CO₂–ATS3 phasing to be different during the two events, with the maximum probabilities reversed in directionality (i.e. with temperature leading at T1 and CO₂ leading at the ACR end, though zero phasing is within 95 % error). Though CH₄ appears to change alongside CO₂ during both intervals, the phasing between CO₂ 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 CO₂ and Antarctic temperature. Bauska et al. (2016), for example, hypothesize that an earlier rise of CO₂ 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 CO₂ 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). CO₂ 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 CO₂ 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 CO₂ series. The addition of the WD paleotemperature record and removal of the Vostok record from ATS3, the updated atmospheric CO₂ 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 CO₂ and Antarctic temperature during the last deglacial warming. We refine the results of these studies by using the high-resolution CO₂ 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 CO₂ and Antarctic climate is small but variable, with phasing ranging from a centennial-scale CO₂ 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 CO₂ 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 CO₂ 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 CO₂ 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 CO₂ measurements through the beginning of the Holocene. Finally, the continued measurement of high-resolution ice core CO₂ records is essential to understand the relationship between CO₂ and global and regional temperature during the last 800 000 years.