Deglacial carbon cycle changes observed in a compilation of 127 benthic δ¹³C time series (20–6 ka)

We present a compilation of 127 time series δ¹³C records from Cibicides wuellerstorfi spanning the last deglaciation (20–6 ka) which is well-suited for reconstructing large-scale carbon cycle changes, especially for comparison with isotope-enabled carbon cycle models. The age models for the δ¹³C records are derived from regional planktic radiocarbon compilations (Stern and Lisiecki, 2014). The δ¹³C records were stacked in nine different regions and then combined using volume-weighted averages to create intermediate, deep, and global δ¹³C stacks. These benthic δ¹³C stacks are used to reconstruct changes in the size of the terrestrial biosphere and deep ocean carbon storage. The timing of change in global mean δ¹³C is interpreted to indicate terrestrial biosphere expansion from 19–6 ka. The δ¹³C gradient between the intermediate and deep ocean, which we interpret as a proxy for deep ocean carbon storage, matches the pattern of atmospheric CO₂ change observed in ice core records. The presence of signals associated with the terrestrial biosphere and atmospheric CO₂ indicates that the compiled δ¹³C records have sufficient spatial coverage and time resolution to accurately reconstruct large-scale carbon cycle changes during the glacial termination.

1 Introduction

On glacial–interglacial timescales, carbon cycle changes redistribute the amount of carbon stored in the deep ocean, atmosphere, and terrestrial biosphere (e.g., Broecker, 1982; Siegenthaler et al., 2005). For example, as atmospheric CO₂ increased across the deglaciation, atmospheric δ¹³C decreased, likely due to the ventilation of respired ¹³C-depleted carbon from the deep ocean (e.g., Schmitt et al., 2012; Eggleston et al., 2016). However, identifying the biogeochemical mechanisms associated with these carbon transfers is complicated by a variety of carbon cycle feedbacks (e.g., Archer et al., 2000; Sigman and Boyle, 2000; Peacock et al., 2006; Toggweiler et al., 2006; Kohfeld and Ridgwell, 2009; Brovkin et al., 2012; Menviel et al., 2012; Galbraith and Jaccard, 2015; Buchanan et al., 2016). This study seeks to improve our understanding of glacial–interglacial carbon cycle changes by reconstructing changes in mean ocean δ¹³C and its vertical gradient and comparing the results with changes in the terrestrial biosphere and atmospheric CO₂.

The δ¹³C of benthic foraminiferal calcite is a well-established carbon cycle proxy, which records the δ¹³C signature of dissolved inorganic carbon (DIC) in seawater at seafloor depths (e.g., Woodruff and Savin, 1985; Zahn et al., 1986; Lutze and Thiel, 1989; Duplessy et al., 1988; Mackensen, 2008; Gottschalk et al., 2016; Schmittner et al., 2017). Averages of benthic foraminiferal δ¹³C time series, called stacks, can improve the signal-to-noise ratio of regional or global seawater changes (e.g., Lisiecki et al., 2008; Lisiecki, 2014). Global mean benthic δ¹³C change is likely caused by changes in terrestrial organic carbon storage (Shackleton, 1977; Curry et al., 1988; Duplessy et al., 1988; Ciais et al., 2012; Peterson et al., 2014), while vertical δ¹³C gradients may record changes in deep ocean carbon storage and atmospheric CO₂ (Oppo and Fairbanks, 1990; Flower et al., 2000; Hodell et al., 2003b; Lisiecki, 2010). The vertical δ¹³C gradient between the surface (high δ¹³C) and deep ocean (low δ¹³C) primarily results from the accumulation of low-δ¹³C-respired organic carbon in deep water, which temporarily sequesters it from the atmosphere. Conversely, vertical mixing of the ocean will tend to ventilate deep ocean carbon to the surface ocean and atmosphere while simultaneously decreasing the vertical δ¹³C gradient. Therefore, the vertical δ¹³C gradient likely records changes in deep ocean carbon storage, which is an important factor controlling glacial–interglacial changes in atmospheric CO₂ (e.g., Schmitt et al., 2012; Eggleston et al., 2016).

Here we compile and analyze 127 high-resolution benthic δ¹³C records from the Atlantic, Pacific, and Indian oceans spanning the last deglaciation to investigate changes in both the ocean and terrestrial biosphere components of the global carbon cycle. Benthic δ¹³C records are combined into regional stacks, which are then used to construct time series of volume-weighted global mean δ¹³C and the vertical δ¹³C gradient between intermediate and deep waters.

We analyze these stacks to test the following hypotheses:

1. The deglacial pattern of global mean ocean δ¹³C change is a proxy for changes in the size of the terrestrial biosphere. If so, global mean δ¹³C should continue to increase after atmospheric CO₂ levels plateau at 11 ka due to the slower response times for ice sheet retreat and ecosystem change (e.g., Hoogakker et al., 2016; Davies-Barnard et al., 2017). We compare the reconstructed global mean δ¹³C change with several carbon cycle model estimates of terrestrial biosphere change. Additionally, we evaluate whether deep Pacific δ¹³C correlates with global mean δ¹³C change as previously assumed (Shackleton et al., 1983; Curry and Oppo, 1997; Lisiecki et al., 2008). This study provides the first opportunity to compare time series of deep Pacific δ¹³C with a volume-weighted global mean δ¹³C stack.

2. Changes in the vertical δ¹³C gradient should closely resemble time series of atmospheric CO₂ if the deglacial CO₂ increase is caused by a decrease in deep ocean carbon storage. This hypothesis is supported by findings on orbital timescales using a smaller number of sites (Oppo and Fairbanks, 1990; Flower et al., 2000; Hodell et al., 2003b; Lisiecki, 2010), but the link between the vertical δ¹³C gradient and CO₂ has not yet been evaluated at millennial timescales or using a global data compilation. Observing such a link would improve our understanding of deglacial atmospheric CO₂ increase and, furthermore, demonstrate that the data compilation presented here has adequate spatial and temporal resolution with sufficiently precise age models to reconstruct millennial-scale changes in the global carbon cycle.

2 Background

2.1 Benthic δ¹³C reconstructions

Measurements of δ¹³C from the calcite tests of epibenthic foraminifera Cibicides wuellerstorfi and related species (Schweizer et al., 2009) are commonly used to trace the spatial distribution of nutrients and deep water masses as well as changes in ocean carbon cycling (e.g., Curry et al., 1988; Duplessy et al., 1988; Curry and Oppo, 2005; Schmittner et al., 2017). Benthic δ¹³C is also slightly influenced (<15 %) by changes in carbonate ion concentration of sea water (Schmittner et al., 2017). Additionally, the Cibicides species C. kullenbergi and C. mundulus, often measured in deep South Atlantic cores, appear to record more depleted δ¹³C values than C. wuellerstorfi (Gottschalk et al., 2016).

Mean δ¹³C has been estimated for the Last Glacial Maximum (LGM, 20 ka) and Late Holocene (6–0 ka) using global compilations of Cibicides wuellerstorfi δ¹³C records (e.g., Shackleton, 1977; Duplessy et al., 1988; Curry et al., 1988; Boyle, 1992; Matsumoto and Lynch-Stieglitz, 1999; Curry and Oppo, 2005; Herguera et al., 2010; Oliver et al., 2010; Hesse et al., 2011; Peterson et al., 2014; Gebbie et al., 2015). These time slice studies include as many as 500 core sites, but generally undersample portions of the ocean with poor carbonate preservation, low primary productivity, and low sedimentation rates (i.e., the Southern Ocean south of 55^∘ S, the Indian Ocean, and the Pacific Ocean). In contrast, some portions of the Atlantic, especially the North Atlantic, are relatively well-sampled with abundant, well-preserved C. wuellerstorfi. Therefore, whole-ocean mean δ¹³C change is less well-constrained than Atlantic δ¹³C.

Because deglacial carbon cycle changes occurred on millennial to centennial timescales (Marcott et al., 2014), observing these changes in the ocean requires a global compilation of high-resolution benthic δ¹³C time series on a consistent age model across the glacial termination. Previous global compilations of δ¹³C time series focus on orbital-scale responses because their age models are not precise enough to analyze the relative timing of carbon cycle changes during the deglaciation (e.g., Lisiecki et al., 2008). For example, Oliver et al. (2010) caution that their global δ¹³C data synthesis, which includes 258 records from many benthic and planktic foraminifera species, should not be used to analyze δ¹³C changes on timescales of less than 10 kyr due to age model uncertainty and the inclusion of low-resolution records. Instead, studies of δ¹³C change across the last glacial termination often use local or regional depth transects that contain high-resolution δ¹³C records with good age control (e.g., Sarnthein et al., 1994; Thornalley et al., 2010; Hoffman and Lund, 2012; Tessin and Lund, 2013; Lund et al., 2015; Oppo et al., 2015; Sikes et al., 2016). In modeling studies, transient simulations are typically compared to a small number of individual benthic δ¹³C records or regional syntheses, presumably due to the limitations of available global δ¹³C compilations (e.g., Köhler et al., 2005, 2010; Brovkin et al., 2007).

2.2 Terrestrial biosphere and mean ocean δ¹³C

A portion of the additional carbon released from the deep ocean since the LGM was taken up by the terrestrial biosphere. The transfer of carbon between the terrestrial biosphere and the deep ocean affects the global mean value of benthic δ¹³C because the mean δ¹³C signature of the terrestrial biosphere is significantly more negative (approximately −25 ‰) than mean ocean δ¹³C (approximately 0 ‰) (Shackleton, 1977). The change in global mean benthic δ¹³C between the LGM and the Holocene is estimated to be 0.32 ‰ ± 0.20 ‰ (Peterson et al., 2014; Gebbie et al., 2015), but the timing of mean benthic δ¹³C change across the deglaciation is not well known.

The amount of terrestrial carbon storage change (soils and vegetation) can be reconstructed in many ways, including terrestrial vegetation proxies and archives (e.g., pollen, paleovegetation); carbon cycle models (e.g., box models, inverse methods, dynamic global vegetation models, biomization methods, etc.); and proxies such as benthic δ¹³C, triple oxygen isotopes (Landais et al., 2007), and atmospheric carbonyl sulfide (Aydin et al., 2016). These methods produce estimates of change in terrestrial carbon storage between the LGM and Holocene varying from 200–1900 PgC due to uncertainties and assumptions associated with each method (see discussion and citations within Peterson et al., 2014).

Due to uncertainties in the total magnitude of change, here we focus on comparing the timing of changes in terrestrial carbon storage and global mean benthic δ¹³C. Models simulate rapid increases in terrestrial carbon storage from approximately 19–10 ka, followed by more gradual changes from 10–0 ka (Kaplan et al., 2002; Joos et al., 2004; Köhler et al., 2005). More recently, the potential effects of changes in poorly constrained carbon reservoirs (e.g., beneath ice sheets and on continental shelves) were evaluated using deglacial simulations of biogeophysical and land carbon changes from the HadCM3 general circulation model (GCM). The model simulated a rapid increase in terrestrial carbon storage from 20–14 ka, different responses between 14–11 ka depending on the model scenario, and then steady, gradual change from 11–4 ka (Davies-Barnard et al., 2017).

Estimates of global mean benthic δ¹³C are also used to remove global changes from individual δ¹³C records to identify patterns of local or regional change, e.g., related to ocean circulation. Because estimates of global mean δ¹³C have only been available for the LGM and Holocene, some studies use deep Pacific δ¹³C time series as a proxy for global mean δ¹³C change (Shackleton et al., 1983; Curry and Oppo, 1997; Lisiecki et al., 2008). Given the large volume and carbon storage capacity of the deep Pacific, its δ¹³C change should be similar in magnitude and timing to the mean ocean δ¹³C change; however, no study has yet confirmed this relationship. For example, low sedimentation rates and poor carbonate preservation in the deep Pacific may limit how well deep Pacific δ¹³C time series resolve changes in mean ocean δ¹³C. Additionally, large changes in Atlantic or Indian Ocean δ¹³C could alter the timing of global mean δ¹³C relative to the Pacific. By constructing a global benthic δ¹³C stack, we can now directly compare deep Pacific δ¹³C with global mean δ¹³C change across the deglaciation.

2.3 Vertical gradients in benthic δ¹³C

A vertical gradient in the δ¹³C of DIC between surface-to-intermediate waters and deep water results from a combination of physical, chemical, and biological processes. The air–sea gas exchange of CO₂ between the atmosphere and surface ocean generates a temperature-dependent fractionation (Lynch-Stieglitz et al., 1995). Biological productivity in the surface ocean preferentially incorporates ¹²C into organic molecules, leaving ¹³C-enriched DIC in surface waters. Conversely, deep water becomes depleted in ¹³C due to remineralization of sinking organic carbon with a δ¹³C signature of approximately −25 ‰. The accumulation of respired organic carbon in the deep ocean gradually increases deep water's DIC concentration while decreasing its δ¹³C value. Thus, sinking organic carbon simultaneously creates vertical gradients in both δ¹³C and DIC, creating low δ¹³C and high DIC in the deep ocean and high δ¹³C and low DIC in the surface ocean. However, deep water δ¹³C is also affected by the transport of relatively high-δ¹³C North Atlantic Deep Water into the deep Atlantic, where it mixes with low-δ¹³C waters from the Southern Ocean (Talley, 2013).

Numerous δ¹³C records from the well-characterized Atlantic Ocean demonstrate an enhanced vertical δ¹³C gradient between intermediate and deep water during the LGM (e.g., Curry and Lohmann, 1982a; Curry et al., 1988; Duplessy et al., 1988; Sarnthein et al., 1994; Hodell et al., 2003b; Curry and Oppo, 2005; Marchitto and Broecker, 2006; Herguera et al., 2010). The less well-sampled Pacific and Indian oceans also show signs of enhanced stratification at the LGM based on stronger vertical δ¹³C gradients and other nutrient and ventilation proxies (e.g., Kallel et al., 1988; Matsumoto and Lynch-Stieglitz, 1999; Matsumoto et al., 2002; Herguera et al., 2010; Lund et al., 2011b; Allen et al., 2015; Sikes et al., 2016).

Figure 1 Locations of 127 core sites compiled for this study, color coded by LGM δ¹³C estimates at each core site. Markers indicate locations of cores in the nine regions: INA is the intermediate North Atlantic; UDNA is the upper deep North Atlantic; LDNA is the lower deep North Atlantic; ISA is the intermediate South Atlantic; DSA is the deep South Atlantic; II is the intermediate Indian; DI is the deep Indian; IP is the intermediate Pacific; DP is the deep Pacific.

Multiple causes have been proposed for stronger vertical δ¹³C gradients during the LGM, including increased surface productivity and export, increased ocean stratification, and changes in preformed δ¹³C in regions of deep water formation (e.g., Matsumoto et al., 2002; Curry and Oppo, 2005; Marchitto and Broecker, 2006; Lynch-Stieglitz et al., 2007; Marinov et al., 2008a, b; Herguera et al., 2010; Hesse et al., 2011; Lund et al., 2011a, b; Allen et al., 2015; Gebbie et al., 2015; Schmittner and Somes, 2016; Gloege et al., 2017; Menviel et al., 2017). Therefore, the large vertical δ¹³C gradient at the LGM could indicate a strong biological pump and/or weak vertical mixing, either of which would increase deep ocean carbon storage. Although studies do not agree about the relative importance of different mechanisms in creating this vertical gradient, the consensus is that the enhanced vertical δ¹³C gradient at the LGM is consistent with greater deep ocean carbon storage and that this carbon was transferred to the atmosphere and terrestrial biosphere during the glacial termination. Multiple processes likely contribute to the deglacial pCO₂ rise (Bauska et al., 2016), including ocean temperature increase, enhanced Southern Ocean mixing rates and the role of sea ice (e.g., Fraçnois et al., 1997; Crosta and Shemesh, 2002; Gildor et al., 2002; Hodell et al., 2003b; Paillard and Parrenin, 2004), decreased alkalinity and carbon inventories (Yu et al., 2014; Kerr et al., 2017), reduced biological pump (Buchanan et al., 2016), enhanced global ocean circulation (Buchanan et al., 2016), and coral reef growth (e.g., Vecsei and Berger, 2004).

On orbital timescales, changes in the intermediate-to-deep vertical δ¹³C gradient closely match atmospheric CO₂, with weaker vertical δ¹³C gradients corresponding to higher CO₂ levels (Oppo and Fairbanks, 1990; Flower et al., 2000; Hodell et al., 2003b; Köhler et al., 2010; Lisiecki, 2010). This relationship suggests that many of the processes affecting CO₂ also alter the vertical δ¹³C gradient. Here we evaluate the relationship between atmospheric CO₂ and vertical δ¹³C change at millennial resolution across the deglaciation. It is beyond the scope of this study to evaluate how much of the change in CO₂ and the vertical δ¹³C gradient at the LGM is associated with specific processes, such as changes in the biological pump (Archer et al., 2003; Köhler et al., 2005; Brovkin et al., 2007; Galbraith and Jaccard, 2015), deep water formation (McManus et al., 2004; Curry and Oppo, 2005), and/or Southern Ocean stratification (Lund et al., 2011b; Burke and Robinson, 2012).

3 Data

This study presents a compilation of 127 previously published benthic δ¹³C time series of Cibicides wuellerstorfi in per mil relative to Vienna PeeDee Belemnite (V.P.D.B.; Fig. 1; Table A1 in the Appendix). Each record in the compilation spans the time range 20–6 ka. Analysis does not extend after 6 ka because cores from several data-sparse regions were either of a too low resolution or missing sediment from 6–0 ka. We only include δ¹³C records with mean sample spacing better than 3 kyr, and 87 % have a mean sample spacing of less than 2 kyr. We excluded any records with sample gaps of 4 kyr or larger and excluded any cores affected by the phytodetritus effect (“Mackensen effect”) as assessed by the original authors and the criteria from Peterson et al. (2014). We included one C. kullenbergi record from the deep South Atlantic (MD07-3076Q; Waelbroeck et al., 2011), which may record a more negative δ¹³C value than C. wuellerstorfi at the LGM (Gottschalk et al., 2016). Additionally, we use some cores with samples labeled “C. spp” that may include some C. kullenbergi (Table A1).

4 Methods

4.1 Age models

For nearly all cores we use the age models of Stern and Lisiecki (2014) based on regional benthic δ¹⁸O alignments and seven regional age models. Each of the seven regions has an age model based on planktic ¹⁴C measurements from multiple cores; ¹⁴C dates are combined across cores by assuming that benthic δ¹⁸O is synchronous within each region (but not necessarily between regions). The first step of this process was generating an initial radiocarbon age model for each of the 61 cores by using that core's radiocarbon dates, the Bayesian age modeling software Bacon (Blaauw and Christen, 2011), the Marine13 calibration (Reimer et al., 2013), and constant 405 ¹⁴C-yr reservoir ages. Bacon was used to estimate ¹⁴C-based ages at specified depths throughout each core, including Monte Carlo uncertainty estimates that increase with distance from the ¹⁴C measurements. To identify the core-specific depths for which ¹⁴C-based ages would be combined, each core's benthic δ¹⁸O record was aligned to an Atlantic or Pacific target core using the alignment software Match (Lisiecki and Lisiecki, 2002). Creating regional age models maximizes the total number of ¹⁴C dates which contribute to each age model. For example, the intermediate Pacific age model is derived from 14 sediment cores that include a total of 160 radiocarbon dates. The final age model for each core in Stern and Lisiecki (2014) was produced by converting from a (transitional) target age model based on benthic δ¹⁸O alignment to a regional composite radiocarbon age model.

Our compilation also includes δ¹³C records from 10 South Atlantic cores that were not included in Stern and Lisiecki (2014) and for which we used the core's published radiocarbon age models (Sortor and Lund, 2011; Hoffman and Lund, 2012; Tessin and Lund, 2013; Lund et al., 2015). These cores are denoted with asterisks in Table A1. The bulk of data compilation work for this study occurred in 2010–2015, and more recently published data are not included.

Figure 2 The three-dimensional structure of δ¹³C in the Indian and Pacific oceans (a, c) and Atlantic (b, d) shown as zonally collapsed cross sections (latitude vs. modern water depth) with the same marker scheme as in Fig. 1. Dotted lines indicate region boundaries. Colors show the δ¹³C value at each site for the LGM (20 ka, a, b) and Holocene (6 ka, c, d). Note that latitudes on the x axis are oriented so the Southern Ocean is in the center of the figure. Additional time slices (in 1 kyr intervals from 20–6 ka) and an animation of deglacial δ¹³C changes can be found in the Supplement.

Table 1 Regional stack information. The total volume represented by the global stack of all nine regions (spanning 0.5–5 km and excluding shallow inland seas, the Southern Ocean, and the Arctic Ocean) is 77.7 % of the whole ocean. The volume of each region is listed as a percent of the global stack volume (rather than whole ocean volume). For each stack we list its δ¹³C value at the LGM (20 ka), Holocene (6 ka), and the Holocene δ¹³C minus LGM δ¹³C difference. The 95 % confidence interval for the global mean δ¹³C change is provided in parentheses. The full time series for each stack is provided in the Supplement.

^a Volume of all regions as a proportion of whole-ocean volume. ^b Atlantic and Pacific Ocean regions, excluding the Indian Ocean regions. ^c Volume-weighted δ¹³C values. ^d Atlantic, Indian, and Pacific Ocean regions.

Stern and Lisiecki (2014) estimate 95 % confidence intervals for each regional age model using 10 000 Monte Carlo age samples for each core from Bacon. Age uncertainty estimates for each region include the effects of any errors in benthic δ¹⁸O alignment because alignment errors would increase scatter in the compiled radiocarbon dates (by aligning portions of cores with different ages) and, thus, increase the observed spread in age estimates. For the time range of 6–20 kyr used in our δ¹³C compilation, the 95 % confidence interval widths of the regional age models range from 0.5–2.0 kyr. Although Match does not quantify alignment uncertainty, alignment uncertainties have been estimated using a similar algorithm, called HMM-Match (Lin et al., 2014). For age models generated either by δ¹⁸O alignment or radiocarbon dating, the amount of age uncertainty depends on the time resolution of the δ¹⁸O or ¹⁴C data, respectively. A comparison of 15 low-latitude Pacific cores found that ¹⁴C-based age uncertainty is comparable to, if not greater than, the uncertainty associated with δ¹⁸O alignments by HMM-Match (Khider et al., 2017).

Figure 3 Regional stacks for the (a) Atlantic and (b) Indian and Pacific oceans. Note the x and y axes are identically scaled. INA is the intermediate North Atlantic; UDNA is the upper deep North Atlantic; LDNA is the lower deep North Atlantic; ISA is the intermediate South Atlantic; DSA is the deep South Atlantic; II is the intermediate Indian; DI is the deep Indian; IP is the intermediate Pacific; DP is the deep Pacific.

4.2 Stacking

After compiling all 127 records of their previously published age models, we use spatial patterns in benthic δ¹³C to define nine ocean regions, for example, based on different LGM δ¹³C values for intermediate and deep sites (Fig. 2). In the North Atlantic, we separate the intermediate North Atlantic (INA, 0.5–2 km) from the upper deep North Atlantic (UDNA, 2–4 km) and the lower deep North Atlantic (LDNA, >4 km). Because the South Atlantic has fewer records than the North Atlantic (Table 1) and a different vertical δ¹³C structure (Fig. 2), we define the intermediate South Atlantic (ISA) as 0.5–2.5 km, and the deep South Atlantic (DSA) as >2.5 km. Although a zonal gradient is evident in the intermediate South Atlantic (Fig. 1), as also observed by Peterson et al. (2014), we combine all ISA records into a single region because only three sites are available in the east. We separate the Indo-Pacific into four regions: the intermediate Indian (II, 0.5–2 km), intermediate Pacific (IP, 0.5–2 km), deep Indian (DI, >2 km), and deep Pacific (DP, >2 km). The longitude boundaries between the Atlantic, Indian, and Pacific basins are the same as in Peterson et al. (2014). Most regions contain at least six sites; however, the intermediate and deep Indian regions each contain only two sites.

Although sea level rises globally by about 130–134 m (Lambeck et al., 2014; Clark et al., 2009) across the deglaciation, the volumetric change associated with deglacial sea level rise is small, less than 3 %. Therefore, we use modern water depths and volumes for all sites at all time steps. This preserves the spatial dimensions of the regions and prevents cores near region boundaries from switching between regions during the deglaciation.

To create regional stacks, we interpolated all benthic δ¹³C records to an even 1 kyr spacing and averaged all records within each region (Fig. 2; Table A1). The δ¹³C time series for sites KNR159-5-90GGC and KNR159-5-22GGC do not include 20 and 6 ka, respectively (Lund et al., 2015); at these times, the relevant sites were excluded from the regional average (Fig. 2 and the animation in the Supplement). Intermediate, deep, and global mean δ¹³C stacks (Figs. 3, 4) are calculated by averaging the regional stacks using volume weighting as a percent of total volume over a depth range of 0.5–5 km (Table 1). Thus, we represent regions proportional to their volume rather than over-representing well-sampled regions.

Figure 4 Volume-weighted stacks. (a) The volume-weighted global stack is calculated based on all nine regional stacks. However, the intermediate and deep stacks shown here only include Atlantic and Pacific data due to the small amount of Indian data. A comparison of these stacks with ones which include the Indian stacks is provided in Fig. S1. (b) Stack uncertainty as characterized by the 95 % confidence interval half-width of each stack.

4.3 Stack limitations and uncertainty

Although our global stack includes benthic δ¹³C records from the Atlantic, Indian, and Pacific oceans, it does not include data from the Southern Ocean, Arctic Ocean, or shallow inland seas. Additionally, our compilation only includes benthic δ¹³C records from below 0.5 km. Although planktic δ¹³C data suggest that mixed layer δ¹³C values may closely track atmospheric δ¹³C change (Eggleston et al., 2016; Hertzberg et al., 2016), we refrain from interpreting or making assumptions about δ¹³C above 0.5 km. It is beyond the scope of the current study to quantify stack uncertainty associated with portions of the ocean which lack C. wuellerstorfi δ¹³C time series.

Uncertainty estimates for δ¹³C from C. wuellerstorfi range from 0.1 ‰ (Marchal and Curry, 2008) to 0.22 ‰ (experiments “LW” and “CW” in Schmittner et al., 2017). Accounting for the carbonate ion concentration of seawater can improve the accuracy of benthic δ¹³C (Schmittner et al., 2017), but estimates of carbonate ion changes throughout the deglaciation are scarce. In the absence of carbonate ion data, a linear regression can be used to convert between C. wuellerstorfi δ¹³C and DIC δ¹³C (regressions LW6 and CW6 in Schmittner et al., 2017). However, because our study focuses on the timing of δ¹³C change rather than its amplitude, we present all δ¹³C data using the values originally measured in foraminiferal calcite.

Interpolating the δ¹³C records to an even 1 kyr spacing introduces an additional source of uncertainty in the data. Although combining information from multiple records inherently risks distorting the true ocean state, this risk is counterbalanced by the potential for improved signal-to-noise ratio when estimating regional and global signals. In the Supplement, we provide the original, uninterpolated records for all 127 sites, which could be used for comparison with transient deglacial ocean circulation experiments. Because age model uncertainties are approximately 1–2 kyr (Stern and Lisiecki, 2014) and some of the δ¹³C records analyzed have sample spacings of 2–3 kyr, our interpretation focuses on δ¹³C features with timescales of about 2 kyr or greater. For example, we do not expect to reconstruct abrupt changes associated with the onset of the Bølling-Allerød or with centennial-scale CO₂ change (Marcott et al., 2014).

We estimate stack uncertainty using Monte Carlo simulations that account for the effects of measurement uncertainty and intra-region δ¹³C variability. Specifically, we generate nominal 95 % confidence intervals for the stacks using 10 000 bootstrapped iterations that randomly resample δ¹³C records from each region. During the resampling process, we also simulate δ¹³C measurement uncertainty in each record by adding Gaussian white noise with a standard deviation of 0.20 ‰ (Gebbie et al., 2015). Multiple runs of our Monte Carlo simulations, each with 10 000 iterations, produce differences in the global benthic δ¹³C stack on the order of 0.02 ‰ at the LGM (20–19 ka) and less during the Holocene.

4.4 Comparison to atmospheric CO₂

To compare the δ¹³C data to atmospheric CO₂ changes from 20–6 ka, we calculate a vertical δ¹³C gradient (Δδ¹³C_I−D) as the difference between the volume-weighted intermediate and deep regional stacks from the Atlantic and Pacific. The Indian Ocean regional stacks are excluded from this vertical gradient calculation because each Indian region includes only two sites, making the Indian regional stacks more susceptible to noise. A global vertical δ¹³C gradient that includes the Atlantic, Indian, and Pacific oceans (AIP Δδ¹³C_I−D) is provided in the Supplement and in the Appendix (Fig. S1, Table A2). Additionally, we evaluate an alternate gradient, $\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{(\mathrm{INA}/\mathrm{2})-\mathrm{DP}}$, defined as the difference between half the intermediate North Atlantic stack and the deep Pacific stack; Lisiecki (2010) found that this gradient optimized correlation with CO₂ from 0–800 ka.

We interpolate a composite ice core CO₂ record (Köhler et al., 2017) to the same 1 kyr resolution as our benthic δ¹³C stacks and calculate correlation coefficients between CO₂ and the vertical gradient of δ¹³C. Additionally, we examine the potential for differences in the timing of CO₂ and δ¹³C change that could be caused by lags in the climate system or age model uncertainty. We evaluate different potential lags by interpolating the CO₂ record with different time offsets, ranging from +1000 to −1000 years in increments of 100 years. For example, a 100-year lag in CO₂ relative to the vertical δ¹³C gradient would be represented by comparing δ¹³C values at 6, 7, …20 ka with CO₂ values at 5.9, 6.9, …, and 19.9 ka. Conversely, a CO₂ lead of 100 years would be suggested if the correlation between the two is maximized for CO₂ values at 6.1, 7.1, …, 20.1 ka.

Testing the significance of correlations between δ¹³C and CO₂ is complicated by the fact that both time series are autocorrelated, i.e., each data point is highly correlated with the value immediately before or after. To reduce the impact of autocorrelation, we pre-whiten the data by taking the difference between successive 1 kyr samples before calculating the linear correlation and its statistical significance. Our final assessment of the statistical significance of the correlations accounts for the reduction in the number of degrees of freedom in the data associated with pre-whitening and allowing time lags between δ¹³C and CO₂ observations.

5 Results

5.1 Comparison to LGM and Holocene reconstructions

Although our compilation of δ¹³C time series includes fewer core sites than some previous studies of LGM δ¹³C, it preserves the large-scale features of these glacial reconstructions, such as enhanced vertical and meridional Atlantic δ¹³C gradients (Fig. 2; e.g., Curry and Oppo, 2005; Peterson et al., 2014). Vertical δ¹³C gradients at the LGM are strongest in the glacial North Atlantic, closely followed by the glacial South Atlantic (Peterson et al., 2014). The most depleted δ¹³C values in the compilation are from the high-latitude deep South Atlantic during the LGM, possibly due to inclusion of data from C. kullenbergi (Gottschalk et al., 2016). Indo-Pacific δ¹³C values for the LGM are similar to equatorial deep South Atlantic records of the same depth and more depleted than North Atlantic δ¹³C values. However, our compilation lacks Indo-Pacific sites deeper than 3.5 km.

At 6 ka the δ¹³C values in this compilation generally resemble the Holocene compilation of Peterson et al. (2014). Minor differences could result from Peterson et al. (2014) using Holocene data from 6–0 ka and including more sites from the North Pacific sites and from 0.5–1.5 km depth.

5.2 Regional stacks

We create nine regional δ¹³C stacks from 20–6 ka (Fig. 3, Table 1). Six of the regional δ¹³C stacks increase steadily from approximately 20–18 to 6 ka (LDNA, DSA, II, IP, DI, DP). Small deviations in the trends of the Indian δ¹³C stacks are interpreted as noise because these stacks each contain only two δ¹³C records. Three Atlantic regions (INA, ISA, UDNA) show a decrease in δ¹³C from approximately 19–15 ka, followed by an increase from 14–6 ka, as described in previous studies (e.g., Hodell et al., 2008, 2010; Thornalley et al., 2010; Lund et al., 2011a; Tessin and Lund, 2013; Oppo et al., 2015). The UDNA δ¹³C stack has a δ¹³C value between the ISA and LDNA from 20–17 ka, approximately matches the LDNA at 16 ka, and then resembles the ISA stack from 14–6 ka. LDNA δ¹³C is slightly greater than DSA δ¹³C except at 10 ka when the two stacks briefly converge. The DI and DSA δ¹³C values are generally similar across the deglaciation except that the DSA δ¹³C begins increasing at 18 ka while the DI δ¹³C increase begins at 16 ka. The intermediate-depth δ¹³C stacks in the Indian and Pacific oceans are very similar for most of the time interval.

Across the deglaciation, the vertical δ¹³C gradient weakens in the Atlantic, most noticeably in the North Atlantic where the INA-LDNA gradient decreases from 1.20 ‰ at 20 ka to 0.31 ‰ at 6 ka. Vertical gradients in the Indian and Pacific oceans show much less change. The largest spread in δ¹³C values is observed from 20–18 ka, when the intermediate North Atlantic and deep South Atlantic regions differ by 1.50 ‰, a difference which decreases to 0.40 ‰ by 10 ka. The maximum difference between regions at 6 ka is 0.91 ‰ between the intermediate North Atlantic and the deep Indian.

5.3 Volume-weighted stacks and global mean δ¹³C stack

A global mean δ¹³C stack is constructed by volume weighting all nine regional stacks. However, we construct two versions of the intermediate and deep δ¹³C stacks, with and without the Indian stacks, because the Indian regions each contain only two records. Both versions of the intermediate and deep stacks show similar trends, but we focus our analysis on the version that uses only the Atlantic and Pacific regions, which should be less susceptible to noise (Fig. 4, Table 1). Results for the intermediate and deep δ¹³C stacks that include the Indian Ocean are provided in the Supplement and in the Appendix (Fig. S1, Table A2).

The volume-weighted intermediate, deep, and global mean δ¹³C stacks increase across the deglaciation, but the magnitude of change is larger for the deep stack (0.46 ‰) than the intermediate stack (0.24 ‰) (Table 1, Fig. 4). We define the vertical δ¹³C gradient, Δδ¹³C_I−D, as the difference between the volume-weighted intermediate and deep stacks that exclude the data-sparse Indian regions. This gradient has a maximum of 0.41 ‰ at 18 ka and decreases to 0.24 ‰ by 6 ka.

Figure 5 Four HadCM3 simulations of terrestrial carbon storage change (biosphere and soils) to the pre-industrial (blue, Davies-Barnard et al., 2017) compared to our volume-weighted global benthic δ¹³C stack (orange). Global mean δ¹³C change most closely resembles the simulation that releases carbon from under ice sheets to the atmosphere and does not store carbon on continental shelves (solid blue). The two y axes are scaled to illustrate the similarity in the pattern of change across the deglaciation but are not meant to imply that the magnitude of change is equivalent.

The volume-weighted global δ¹³C stack holds nearly steady from 20 to 19 ka at approximately 0.00 ‰ (95 % CI: −0.13 to 0.12 ‰ at 19 ka) and then increases from 18–6 ka, reaching a value of 0.39 ‰ (95 % CI: 0.24 to 0.49 ‰) at 6 ka. The change from 20 to 6 ka in the global stack is 0.36 ‰ (95 % CI: 0.24 to 0.50 ‰), which agrees to within uncertainty with the LGM-to-Holocene δ¹³C change of 0.38 ‰ (95 % CI: 0.30 to 0.46 ‰) estimated for 0.5–5 km by Peterson et al. (2014). Recall that the mean δ¹³C estimate from our global stack is not quite a whole-ocean δ¹³C estimate because we do not include data from the surface (<0.5 km), Southern Ocean (>65^∘ S), or bottom waters (>5 km). Estimates of whole-ocean δ¹³C change are slightly smaller at 0.34 ‰ (95 % CI: 0.15 to 0.53 ‰; Peterson et al., 2014) and 0.32 ‰ (95 % CI: 0.12 to 0.52 ‰; Gebbie et al., 2015) because the surface ocean (0–0.5 km) has less deglacial change (Eggleston et al., 2016; Hertzberg et al., 2016).

6 Discussion

6.1 Terrestrial carbon storage and global mean benthic δ¹³C

The long-standing explanation for mean benthic δ¹³C change across the deglaciation is an increase in the size of the terrestrial biosphere (Shackleton, 1977; Curry et al., 1988; Duplessy et al., 1988). Here we compare the timing of changes in our global mean δ¹³C stack (i.e., a monotonic increase from 19–6 ka, Fig. 4b) with model simulations and other terrestrial biosphere reconstructions.

A carbon isotope-enabled transient model from Lund-Potsdam-Jena Dynamic Global Vegetation Model (LPJ-DGVM) simulated a mean ocean δ¹³C increase beginning at 21 ka, with the most rapid changes occurring from 17–10.5 ka (Kaplan et al., 2002). In these experiments, the terrestrial biosphere began expanding around 18–16.5 ka (Joos et al., 2004; Köhler et al., 2005) and rapidly increased from 17–9 ka, with 70 % of terrestrial carbon storage change occurring before the Holocene (11.5 ka; Kaplan et al., 2002). Similarly, 67 % of the change in our global δ¹³C stack occurs between 19–11 ka while the remaining 33 % occurs from 11–6 ka.

Figure 6 (a) Time series of global mean ocean δ¹³C stack and deep Pacific δ¹³C stacks. (b) Half-width 95 % confidence intervals for the global mean stack and deep Pacific stack. (c) Deep Pacific δ¹³C stack vs. global mean δ¹³C stack. Each point is the δ¹³C value for one time slice with 95 % confidence intervals (vertical and horizontal error bars). Time across the deglaciation progresses toward the upper right corner. The best-fit linear regression is plotted as a solid line with 95 % confidence interval (dashed lines).

Table 2 Correlation coefficients and p values between pre-whitened records. Pre-whitening reduces the impact of autocorrelation in the time series. Calculated p values account for the reduced degrees of freedom in pre-whitened and time-lagged correlations. Non-zero CO₂ time shifts indicate the lead/lag adjustment that maximizes the correlation between atmospheric CO₂ Köhler et al. (2017) and the δ¹³C stack or gradient.

Simulations from HadCM3 estimated that 45 %–70 % of terrestrial biosphere expansion occurred between 18–14 ka (Davies-Barnard et al., 2017). Dramatically different trends were observed from 14–6 ka in simulations with different assumptions about carbon storage under glacial ice sheets and on continental shelves. The simulation that most closely resembles our global mean δ¹³C stack is the simulation that releases carbon from under ice sheets to the atmosphere and does not accumulate carbon on exposed continental shelves (Fig. 5). This simulation is also the only one which agrees with terrestrial carbon storage change estimates of 440 PgC based on whole-ocean mean δ¹³C change (e.g., Peterson et al., 2014).

Holocene simulations using a global pollen synthesis, the biomization method, and models of climate and vegetation (HadCM3, FAMOUS, and BIOME4) suggest that the global average area for most carbon-rich megabiomes (i.e., excluding grasslands and dry shrubland) increased from 10–2 ka and net primary productivity increased from 8–2 ka (Hoogakker et al., 2016). This is consistent with our observation that the global mean benthic δ¹³C trend continued until at least 6 ka. Dramatic land use changes from agricultural practices, another potential mechanism for terrestrial carbon change, did not begin until 4.5 ka (Ruddiman and Ellis, 2009). More detailed evaluation of Holocene terrestrial carbon storage changes will require improved spatial coverage for δ¹³C records from 6–0 ka.

6.2 Deep Pacific and global mean δ¹³C

Previous studies have assumed deep Pacific δ¹³C can be used as a proxy for global mean δ¹³C because the deep Pacific constitutes about 30 % of the ocean volume and is not strongly affected by shifting water mass boundaries (e.g., Shackleton et al., 1983; Curry and Oppo, 1997; Lisiecki et al., 2008). From 20–6 ka, the global mean and deep Pacific δ¹³C stacks show similar patterns of change (Fig. 6) and fall along a tight regression line

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{global}}=\mathrm{1.04}\pm \mathrm{0.06}\mathrm{‰}\times {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DP}}+\mathrm{0.19}\pm \mathrm{0.01}\mathrm{‰}\end{array}$$

The two time series are highly correlated (r²=0.99), which is not surprising because the large volume of the deep Pacific exerts a strong influence on the global mean δ¹³C stack. When the stacks are pre-whitened to account for autocorrelation (Table 2), their correlation is weaker (r²=0.46) but statistically significant (p=0.05).

Figure 7 Comparison of atmospheric CO₂ with (a) the vertical δ¹³C gradient (excluding the data-sparse Indian Ocean regions) and (b) $\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{(\mathrm{INA}/\mathrm{2})-\mathrm{DP}}$. Both gradients have a statistically significant correlation with atmospheric CO₂ records (red circles, Köhler et al., 2017). The y axes for (a) and (b) are scaled differently because $\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{(\mathrm{INA}/\mathrm{2})-\mathrm{DP}}$ scales the intermediate North Atlantic stack by half before subtracting the deep Pacific stack while the δ¹³C gradient is the difference between the intermediate and deep stacks (Atlantic and Pacific oceans only).

Alternatively, a carbon cycle box model simulated a strong correlation between deep Pacific δ¹³C and CO₂ across several glacial cycles (r²=0.96; Köhler et al., 2010). The correlation between CO₂ and our deep Pacific δ¹³C stack is statistically significant after pre-whitening (r²=0.57, p=0.02), but global mean δ¹³C and CO₂ are not (r²=0.28, p=0.16). Our compilation of Pacific records is likely insufficient to determine whether deep Pacific δ¹³C correlates better with global mean δ¹³C or atmospheric CO₂. This issue could be better resolved using a δ¹³C compilation spanning multiple glacial cycles and including more deep Pacific sites.

6.3 Vertical δ¹³C gradient and atmospheric CO₂

The vertical δ¹³C gradient (Δδ¹³C_I−D) in our compilation resembles the inverse of CO₂ change across the deglaciation (Fig. 7), as would be expected if they are both strongly influenced by changes in deep ocean carbon storage (Flower et al., 2000; Oppo and Horowitz, 2000; Hodell et al., 2003b). Alternatively, one orbital-scale study found a stronger correlation with CO₂ using the gradient between the deep Pacific and half the INA δ¹³C stack ($\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{(\mathrm{INA}/\mathrm{2})-\mathrm{DP}}$), Lisiecki (2010). Both vertical δ¹³C gradients (Δδ¹³C_I−D and $\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{(\mathrm{INA}/\mathrm{2})-\mathrm{DP}}$) decrease from 18–11 ka over the same time interval that atmospheric CO₂ increases. In contrast, the global mean δ¹³C stack increases at a relatively steady pace from 19–6 ka. Thus, the δ¹³C gradients record a distinctly different signal than global mean δ¹³C.

The δ¹³C gradients decrease most rapidly across two time steps, 18–17 and 12–11 ka. The first change at 18 ka is approximately synchronous with the start of atmospheric CO₂ rise (Marcott et al., 2014; Köhler et al., 2017) and a decrease of 0.3 ‰ in the δ¹³C of atmospheric CO₂ (Eggleston et al., 2016). In the Southern Ocean at 18 ka, proxy records indicate a decrease in aeolian dust deposition accompanied by lower marine productivity (Martínez-Garcia et al., 2009) and a decrease in winter sea ice cover, which likely reduced vertical stratification (Ferrari et al., 2014). The second rapid change in the vertical δ¹³C gradients at 12 ka approximately coincides with rapidly increasing atmospheric CO₂ from 13–11.5 ka and a decrease of 0.1 ‰ in the δ¹³C of atmospheric CO₂.

From 11 to 6 ka, atmospheric CO₂ remains nearly constant with a small (approximately 10 ppm) decrease from 11–8 ka. The vertical δ¹³C gradients are also relatively steady from 11–6 ka, with a slight increase in both gradients from 9–8 ka. The small decrease in atmospheric CO₂ beginning at 11 ka (Marcott et al., 2014) has been variously attributed to growth of the terrestrial biosphere, sea level rise, and an increase in gas exchange through reduced sea ice cover (Kaplan et al., 2002; Joos et al., 2004; Köhler and Fischer, 2004; Köhler et al., 2005, 2010).

Although CO₂ correlates strongly with both Δδ¹³C_I−D ($r^{\mathrm{2}}=-\mathrm{0.96}$) and $\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{(\mathrm{INA}/\mathrm{2})-\mathrm{DP}}$ ($r^{\mathrm{2}}=-\mathrm{0.98}$), we must pre-whiten these time series to remove autocorrelation before assessing the statistical significance of their correlation. At the 95 % confidence level, atmospheric CO₂ significantly correlates with both Δδ¹³C_I−D ($r^{\mathrm{2}}=-\mathrm{0.51}$; p=0.03) and $\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{(\mathrm{INA}/\mathrm{2})-\mathrm{DP}}$ ($r^{\mathrm{2}}=-\mathrm{0.69}$; p<0.01) (Fig. S1, Table 2). Better correlation with $\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{(\mathrm{INA}/\mathrm{2})-\mathrm{DP}}$ could be because of better age control and higher resolution δ¹³C records in the INA region than the other intermediate regions. Determining whether the $\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{(\mathrm{INA}/\mathrm{2})-\mathrm{DP}}$ gradient or the global vertical δ¹³C gradient correlates better with atmospheric CO₂ will require data with better spatial coverage and/or a longer time span.

Because our comparison of the vertical δ¹³C gradient and CO₂ could be affected by lags within the carbon cycle or age model uncertainty, we additionally investigate whether the correlations between CO₂ and the vertical δ¹³C gradient would be improved by age model shifts (Table 2). The correlation between CO₂ is maximized when Δδ¹³C_I−D or $\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{(\mathrm{INA}/\mathrm{2})-\mathrm{DP}}$ lags CO₂ by 400 years (Table 2), which is within the age uncertainty of the sediment core age models. Thus, changes in atmospheric CO₂ and vertical δ¹³C gradients appear synchronous to within age model uncertainty.

Processes that potentially explain atmospheric CO₂ change during glacial cycles include the efficiency of the biological pump (Martínez-Garcia et al., 2009; Galbraith and Jaccard, 2015), circulation changes (Ferrari et al., 2014; Schmittner and Lund, 2015; Lacerra et al., 2017; Menviel et al., 2017; Sikes et al., 2017; Wagner and Hendy, 2017), or a combination of multiple processes (Bauska et al., 2016; Skinner et al., 2017). Different processes could influence the carbon cycle on different timescales (Bauska et al., 2016; Kohfeld and Chase, 2017) and/or in different regions (e.g., Gu et al., 2017) and complicate interpretations of which processes are most responsible for atmospheric CO₂ change. However, because both productivity and circulation change would affect the vertical δ¹³C gradient while changing atmospheric CO₂, we interpret our results as supporting the importance of the deep ocean as a reservoir for storing glacial carbon related to either or both processes. Furthermore, these results support the use of vertical δ¹³C gradients as a proxy for glacial–interglacial CO₂ change on both orbital and millennial timescales (Lisiecki, 2010).

7 Conclusions

We present regional δ¹³C stacks and volume-weighted intermediate, deep, and global mean δ¹³C stacks from a compilation of 127 benthic C. wuellerstorfi δ¹³C records, which span 20 to 6 kyr with a mean age resolution better than 2 kyr. Age models are based on δ¹⁸O alignments to regional stacks with radiocarbon dating and age model uncertainties of approximately 1–2 kyr. Our compilation shows spatial patterns in benthic δ¹³C that are similar to higher resolution reconstructions of the Holocene and Last Glacial Maximum. The volume-weighted mean δ¹³C change estimated from these 127 records is 0.36 ‰ (95 % CI: 0.24 to 0.50 ‰), similar to the estimate of Peterson et al. (2014) for 0.5–5 km based on 480 records.

Importantly, this global compilation of benthic δ¹³C time series also allows us to evaluate the timing of change in the mean and vertical gradient of δ¹³C and compare them with other carbon cycle changes. The volume-weighted global δ¹³C stack increases from 19 to 6 ka and likely reflects terrestrial biosphere growth, in agreement with model simulations (Kaplan et al., 2002; Joos et al., 2004; Davies-Barnard et al., 2017). To constrain the timing of the end of terrestrial biosphere expansion, future work should focus on extending the global stack through the Late Holocene. Furthermore, δ¹³C changes from 20 to 6 ka suggest that a deep Pacific δ¹³C stack approximates global mean δ¹³C with an offset of 0.19 ‰. Vertical δ¹³C gradients between intermediate and deep water (Δδ¹³C_I−D) and between the intermediate North Atlantic and deep Pacific ($\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{(\mathrm{INA}/\mathrm{2})-\mathrm{DP}}$) are interpreted as proxies for change in deep ocean carbon storage. Millennial-scale features in Δδ¹³C_I−D and $\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{(\mathrm{INA}/\mathrm{2})-\mathrm{DP}}$ are significantly correlated with atmospheric CO₂ changes from 20–6 ka.

Based on these analyses, we conclude that the four-dimensional compilation of globally distributed δ¹³C time series presented here provides useful constraints for global carbon cycle reconstructions and for comparison with deglacial simulations from isotope-enabled Earth system models.

Appendix A

Sarnthein (1994j)Sarnthein (1994g)Sarnthein (2004b)Dorschel et al. (2005)Curry et al. (1999)Venz et al. (1999)Hüls (1999)Freudenthal et al. (2002a)Jung and Sarnthein (2003g)Sarnthein (1994k)Sarnthein (1994i)Raymo et al. (2004)Oppo and Fairbanks (1990)Weinelt and Sarnthein (2003)Demenocal et al. (1992)Curry et al. (1988j)Millo et al. (2006)McIntyre et al. (1999)Raymo et al. (2004)Sarnthein (2004a)Sarnthein (1994c)Duplessy (1997)Oppo et al. (2006)Jung and Sarnthein (2003c)Duplessy (1982)Sarnthein and Tiedemann (1989)Tjallingii (2008)Sarnthein (2004b)Freudenthal et al. (2002b)Mulitza (2009)Jansen and Veum (1990)Sarnthein (1994h)Schönfeld et al. (2003a)Weinelt (2004)Cortijo et al. (1999)Winn and Sarnthein (1991)Sarnthein (1994d)Sarnthein (1994b)Oppo and Lehman (1995)Richter (2001)Sarnthein (1994a)Fronval and Jansen (1997)Jung and Sarnthein (2003b)Curry et al. (1988c)Jung and Sarnthein (2003f)Curry and Lohmann (1983a)Sarnthein (1994f)Curry et al. (1988c)Curry et al. (1988b)Curry et al. (1988m)Skinner and Shackleton (2004)Curry and Lohmann (1983b)Curry et al. (1988h)Curry et al. (1988k)Zarriess and Mackensen (2011)Jung and Sarnthein (2003a)Schönfeld et al. (2003b)Curry and Lohmann (1983f)Curry et al. (1988g)Jung and Sarnthein (2003d)Curry and Lohmann (1983c)Curry et al. (1988i)Curry et al. (1988l)Sarnthein (1994l)Jung and Sarnthein (2003e)Hodell et al. (2008)Labeyrie (1996)Curry et al. (1988d)Curry and Oppo (2005)Curry and Oppo (1997)Curry and Lohmann (1983e)Sarnthein (2004b)Curry et al. (1988f)Sarnthein (1994e)Curry et al. (1988e)Bickert et al. (2003a)Curry and Lohmann (1983d)Curry and Lohmann (1983g)Arz et al. (1999)Curry et al. (1988a)Lund et al. (2015)Lund et al. (2015)Oppo and Horowitz (2000)Lund et al. (2015)Lund et al. (2015)Oppo and Horowitz (2000)Lund et al. (2015)Hodell et al. (2003a)Lund et al. (2015)Curry and Lohmann (1982b)Oppo and Horowitz (2000)Rau et al. (2002)Lund et al. (2015)Sarnthein (1997)Lund et al. (2015)Lund et al. (2015)Schmiedl and Mackensen (1997)Bickert et al. (2003b)Curry (2004)Bickert et al. (2003c)Lund et al. (2015)Hodell et al. (2003b)Oppo and Fairbanks (1987)Waelbroeck et al. (2011)Mackensen et al. (2001b)Lund et al. (2015)Bickert et al. (2001)Bickert et al. (2009a)Bickert et al. (2009b)Hodell et al. (2003b)Bickert et al. (2009a)Holbourn et al. (2005)Schmiedl and Leuschner (2005)Sirocko et al. (2000a)Sirocko et al. (2000b)Bostock et al. (2004)Pahnke and Zahn (2005)Mix et al. (1991)Chen (2003)Wang et al. (1999)Stott et al. (2000)Shackleton et al. (1992)Ohkushi et al. (2003)Tian et al. (2002)Zhang et al. (2007)Imbrie et al. (1992)Mix et al. (1995)Marchitto et al. (2005)

Table A1 Supplemental table of the name, location, region, and reference for each record in this δ¹³C synthesis. Asterisks mark cores for which we use the cited authors' radiocarbon age model instead of using the regional age model from Stern and Lisiecki (2014). Data repository links for each δ¹³C record are provided in Table S1 in the Supplement.

Table A2 Correlation coefficients and p values between records. The upper rows use the raw data, and the bottom rows use pre-whitened data to account for autocorrelated time series. AIP gradients include Atlantic, Indian, and Pacific regions, and AP gradients include only Atlantic and Pacific regions. To investigate possible leads/lags between records, we shift the atmospheric CO₂ record in 100-year increments relative to the δ¹³C stacks and, for brevity, list only the best correlations. All p values account for reduction in degrees of freedom due to pre-whitening and/or time shifting.