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

. We present a compilation of 127 time series δ 13 C records from Cibicides wuellerstorﬁ 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 δ 13 C records are derived from regional planktic radiocarbon compilations (Stern and Lisiecki, 2014). The δ 13 C records were stacked in nine different regions and then combined using volume-weighted averages to create intermediate, deep, and global δ 13 C stacks. These benthic δ 13 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 δ 13 C is interpreted to indicate terrestrial biosphere expansion from 19–6 ka. The δ 13 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 2 change observed in ice core records. The presence of signals associated with the terrestrial biosphere and atmospheric CO 2 indicates that the compiled δ 13 C records have sufﬁcient spatial coverage and time resolution to accurately reconstruct large-scale carbon cycle changes during the glacial termination.


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 2 increased across the deglaciation, atmospheric δ 13 C decreased, likely due to the ventilation of respired 13 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 δ 13 C and its vertical gradient and comparing the results with changes in the terrestrial biosphere and atmospheric CO 2 .
The δ 13 C of benthic foraminiferal calcite is a wellestablished carbon cycle proxy, which records the δ 13 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 δ 13 C time series, called stacks, can improve the signal-to-noise ratio of regional or global seawater changes (e.g., Lisiecki et al., 2008;. Global mean benthic δ 13 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 δ 13 C gradients may record changes in deep ocean carbon storage and atmospheric CO 2 Flower et al., 2000;Hodell et al., 2003b;Lisiecki, 2010). The vertical δ 13 C gradient between the surface (high δ 13 C) and deep ocean (low δ 13 C) primarily results from the accumulation of low-δ 13 Crespired organic carbon in deep water, which temporarily sequesters it from the atmosphere. Conversely, vertical mix-ing of the ocean will tend to ventilate deep ocean carbon to the surface ocean and atmosphere while simultaneously decreasing the vertical δ 13 C gradient. Therefore, the vertical δ 13 C gradient likely records changes in deep ocean carbon storage, which is an important factor controlling glacialinterglacial changes in atmospheric CO 2 (e.g., Schmitt et al., 2012;Eggleston et al., 2016).
Here we compile and analyze 127 high-resolution benthic δ 13 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 δ 13 C records are combined into regional stacks, which are then used to construct time series of volume-weighted global mean δ 13 C and the vertical δ 13 C gradient between intermediate and deep waters.
We analyze these stacks to test the following hypotheses: 1. The deglacial pattern of global mean ocean δ 13 C change is a proxy for changes in the size of the terrestrial biosphere. If so, global mean δ 13 C should continue to increase after atmospheric CO 2 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 δ 13 C change with several carbon cycle model estimates of terrestrial biosphere change. Additionally, we evaluate whether deep Pacific δ 13 C correlates with global mean δ 13 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 δ 13 C with a volume-weighted global mean δ 13 C stack.
2. Changes in the vertical δ 13 C gradient should closely resemble time series of atmospheric CO 2 if the deglacial CO 2 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 Flower et al., 2000;Hodell et al., 2003b;Lisiecki, 2010), but the link between the vertical δ 13 C gradient and CO 2 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 2 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.

Benthic δ 13 C reconstructions
Measurements of δ 13 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 δ 13 C is also slightly influenced (< 15 %) by changes in carbonate ion concentration of sea water . Additionally, the Cibicides species C. kullenbergi and C. mundulus, often measured in deep South Atlantic cores, appear to record more depleted δ 13 C values than C. wuellerstorfi (Gottschalk et al., 2016).
Mean δ 13 C has been estimated for the Last Glacial Maximum (LGM, 20 ka) and Late Holocene (6-0 ka) using global compilations of Cibicides wuellerstorfi δ 13 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 δ 13 C change is less well-constrained than Atlantic δ 13 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 δ 13 C time series on a consistent age model across the glacial termination. Previous global compilations of δ 13 C time series focus on orbitalscale 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 δ 13 C data synthesis, which includes 258 records from many benthic and planktic foraminifera species, should not be used to analyze δ 13 C changes on timescales of less than 10 kyr due to age model uncertainty and the inclusion of lowresolution records. Instead, studies of δ 13 C change across the last glacial termination often use local or regional depth transects that contain high-resolution δ 13 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 δ 13 C records or regional syntheses, presumably due to the limitations of available global δ 13 C compilations (e.g., Köhler et al., 2005Köhler et al., , 2010Brovkin et al., 2007).
C. D. Peterson and L. E. Lisiecki: Deglacial benthic carbon isotope stacks 1231 2.2 Terrestrial biosphere and mean ocean δ 13 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 δ 13 C because the mean δ 13 C signature of the terrestrial biosphere is significantly more negative (approximately −25 ‰) than mean ocean δ 13 C (approximately 0 ‰) (Shackleton, 1977). The change in global mean benthic δ 13 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 δ 13 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 δ 13 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 δ 13 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 δ 13 C are also used to remove global changes from individual δ 13 C records to identify patterns of local or regional change, e.g., related to ocean circulation. Because estimates of global mean δ 13 C have only been available for the LGM and Holocene, some studies use deep Pacific δ 13 C time series as a proxy for global mean δ 13 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 δ 13 C change should be similar in magnitude and timing to the mean ocean δ 13 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 δ 13 C time series resolve changes in mean ocean δ 13 C.
Additionally, large changes in Atlantic or Indian Ocean δ 13 C could alter the timing of global mean δ 13 C relative to the Pacific. By constructing a global benthic δ 13 C stack, we can now directly compare deep Pacific δ 13 C with global mean δ 13 C change across the deglaciation.

Vertical gradients in benthic δ 13 C
A vertical gradient in the δ 13 C of DIC between surface-tointermediate waters and deep water results from a combination of physical, chemical, and biological processes. The air-sea gas exchange of CO 2 between the atmosphere and surface ocean generates a temperature-dependent fractionation (Lynch-Stieglitz et al., 1995). Biological productivity in the surface ocean preferentially incorporates 12 C into organic molecules, leaving 13 C-enriched DIC in surface waters. Conversely, deep water becomes depleted in 13 C due to remineralization of sinking organic carbon with a δ 13 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 δ 13 C value. Thus, sinking organic carbon simultaneously creates vertical gradients in both δ 13 C and DIC, creating low δ 13 C and high DIC in the deep ocean and high δ 13 C and low DIC in the surface ocean. However, deep water δ 13 C is also affected by the transport of relatively high-δ 13 C North Atlantic Deep Water into the deep Atlantic, where it mixes with low-δ 13 C waters from the Southern Ocean (Talley, 2013).
Multiple causes have been proposed for stronger vertical δ 13 C gradients during the LGM, including increased surface productivity and export, increased ocean stratification, and changes in preformed δ 13 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 δ 13 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 www.clim-past.net/14/1229/2018/ in creating this vertical gradient, the consensus is that the enhanced vertical δ 13 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 2 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-todeep vertical δ 13 C gradient closely match atmospheric CO 2 , with weaker vertical δ 13 C gradients corresponding to higher CO 2 levels 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 2 also alter the vertical δ 13 C gradient. Here we evaluate the relationship between atmospheric CO 2 and vertical δ 13 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 2 and the vertical δ 13 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 Curry and Oppo, 2005), and/or Southern Ocean stratification (Lund et al., 2011b;Burke and Robinson, 2012).

Data
This study presents a compilation of 127 previously published benthic δ 13 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 δ 13 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 δ 13 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).

Age models
For nearly all cores we use the age models of Stern and Lisiecki (2014) based on regional benthic δ 18 O alignments and seven regional age models. Each of the seven regions has an age model based on planktic 14 C measurements from multiple cores; 14 C dates are combined across cores by assuming that benthic δ 18 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 14 C-yr reservoir ages. Bacon was used to estimate 14 C-based ages at specified depths throughout each core, including Monte Carlo uncertainty estimates that increase with distance from the 14 C measurements. To identify the core-specific depths for which 14 C-based ages would be combined, each core's benthic δ 18 O record was aligned to an Atlantic or Pacific target core using the alignment software Match ). Creating regional age models maximizes the total number of 14 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 δ 18 O alignment to a regional composite radiocarbon age model.
Our compilation also includes δ 13 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. 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 δ 18 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 δ 13 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 δ 18 O alignment or radiocarbon dating, the amount of age uncertainty depends on the time resolution of the δ 18 O or 14 C data, respectively. A comparison of 15 low-latitude Pacific cores found that 14 C-based age uncertainty is comparable to, if not greater than, the uncertainty associated with δ 18 O alignments by HMM-Match (Khider et al., 2017).

Stacking
After compiling all 127 records of their previously published age models, we use spatial patterns in benthic δ 13 C to define nine ocean regions, for example, based on different LGM δ 13 C values for intermediate and deep sites (Fig. 2). In the North Atlantic, we separate the intermediate North At-lantic (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 δ 13 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 δ 13 C records to an even 1 kyr spacing and averaged all records within each region ( Fig. 2; Table A1). The δ 13 C time series for sites KNR159-5-90GGC and KNR159-5-22GGC do not include 20 and 6 ka, respectively ; 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 δ 13 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.

Stack limitations and uncertainty
Although our global stack includes benthic δ 13 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 δ 13 C records from below 0.5 km. Although planktic δ 13 C data suggest that mixed layer δ 13 C values may closely track atmospheric δ 13 C change (Eggleston et al., 2016;Hertzberg et al., 2016), we refrain from interpreting or making assumptions about δ 13 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 δ 13 C time series.
Uncertainty estimates for δ 13 C from C. wuellerstorfi range from 0.1 ‰ (Marchal and Curry, 2008) to 0.22 ‰ (exper- Figure 2. The three-dimensional structure of δ 13 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 δ 13 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 δ 13 C changes can be found in the Supplement. iments "LW" and "CW" in Schmittner et al., 2017). Accounting for the carbonate ion concentration of seawater can improve the accuracy of benthic δ 13 C , 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 δ 13 C and DIC δ 13 C (regressions LW6 and CW6 in Schmittner et al., 2017). However, because our study focuses on the timing of δ 13 C change rather than its amplitude, we present all δ 13 C data using the values originally measured in foraminiferal calcite.
Interpolating the δ 13 C records to an even 1 kyr spacing introduces an additional source of uncertainty in the data. Although combining information from multiple records inher-ently 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  and some of the δ 13 C records analyzed have sample spacings of 2-3 kyr, our interpretation focuses on δ 13 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 2 change (Marcott et al., 2014). 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 δ 13 C value at the LGM (20 ka), Holocene (6 ka), and the Holocene δ 13 C minus LGM δ 13 C difference. The 95 % confidence interval for the global mean δ 13 C change is provided in parentheses. The full time series for each stack is provided in the Supplement.

Region name
Sites in stack %Volume δ 13 C Hol (‰) δ 13 C LGM (‰) δ 13 C Hol−LGM (‰)  We estimate stack uncertainty using Monte Carlo simulations that account for the effects of measurement uncertainty and intra-region δ 13 C variability. Specifically, we generate nominal 95 % confidence intervals for the stacks using 10 000 bootstrapped iterations that randomly resample δ 13 C records from each region. During the resampling process, we also simulate δ 13 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 dif-ferences in the global benthic δ 13 C stack on the order of 0.02 ‰ at the LGM (20-19 ka) and less during the Holocene.

Comparison to atmospheric CO 2
To compare the δ 13 C data to atmospheric CO 2 changes from 20-6 ka, we calculate a vertical δ 13 C gradient ( δ 13 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 in- cludes only two sites, making the Indian regional stacks more susceptible to noise. A global vertical δ 13 C gradient that includes the Atlantic, Indian, and Pacific oceans (AIP δ 13 C I−D ) is provided in the Supplement and in the Appendix (Fig. S1, Table A2). Additionally, we evaluate an alternate gradient, δ 13 C (INA/2)−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 2 from 0-800 ka.
We interpolate a composite ice core CO 2 record (Köhler et al., 2017) to the same 1 kyr resolution as our benthic δ 13 C stacks and calculate correlation coefficients between CO 2 and the vertical gradient of δ 13 C. Additionally, we examine the potential for differences in the timing of CO 2 and δ 13 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 2 record with different time offsets, ranging from +1000 to −1000 years in increments of 100 years. For example, a 100-year lag in CO 2 relative to the vertical δ 13 C gradient would be represented by comparing δ 13 C values at 6, 7, . . . 20 ka with CO 2 values at 5.9, 6.9, . . . , and 19.9 ka. Conversely, a CO 2 lead of 100 years would be suggested if the correlation between the two is maximized for CO 2 values at 6.1, 7.1, . . . , 20.1 ka.
Testing the significance of correlations between δ 13 C and CO 2 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 free-dom in the data associated with pre-whitening and allowing time lags between δ 13 C and CO 2 observations.

Comparison to LGM and Holocene reconstructions
Although our compilation of δ 13 C time series includes fewer core sites than some previous studies of LGM δ 13 C, it preserves the large-scale features of these glacial reconstructions, such as enhanced vertical and meridional Atlantic δ 13 C gradients ( Fig. 2; e.g., Curry and Oppo, 2005;Peterson et al., 2014). Vertical δ 13 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 δ 13 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 δ 13 C values for the LGM are similar to equatorial deep South Atlantic records of the same depth and more depleted than North Atlantic δ 13 C values. However, our compilation lacks Indo-Pacific sites deeper than 3.5 km.
At 6 ka the δ 13 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.

Regional stacks
We create nine regional δ 13 C stacks from 20-6 ka (Fig. 3, Table 1). Six of the regional δ 13 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 δ 13 C stacks are interpreted as noise because these stacks each contain only two δ 13 C records. Three Atlantic regions (INA, ISA, UDNA) show a decrease in δ 13 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;Hodell et al., 2010;Thornalley et al., 2010;Lund et al., 2011a;Tessin and Lund, 2013;Oppo et al., 2015). The UDNA δ 13 C stack has a δ 13 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 δ 13 C is slightly greater than DSA δ 13 C except at 10 ka when the two stacks briefly converge. The DI and DSA δ 13 C values are generally similar across the deglaciation except that the DSA δ 13 C begins increasing at 18 ka while the DI δ 13 C increase begins at 16 ka. The intermediate-depth δ 13 C stacks in the Indian and Pacific oceans are very similar for most of the time interval.
Across the deglaciation, the vertical δ 13 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 δ 13 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.

Volume-weighted stacks and global mean δ 13 C stack
A global mean δ 13 C stack is constructed by volume weighting all nine regional stacks. However, we construct two versions of the intermediate and deep δ 13 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 δ 13 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 δ 13 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 δ 13 C gradient, δ 13 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.

Terrestrial carbon storage and global mean benthic δ 13 C
The long-standing explanation for mean benthic δ 13 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 δ 13 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 δ 13 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 δ 13 C stack occurs between 19-11 ka while the remaining 33 % occurs from 11-6 ka.
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 δ 13 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 δ 13 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-  Davies-Barnard et al., 2017) compared to our volume-weighted global benthic δ 13 C stack (orange). Global mean δ 13 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. 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 2 time shifts indicate the lead/lag adjustment that maximizes the correlation between atmospheric CO 2 Köhler et al. (2017)  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 δ 13 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 δ 13 C records from 6-0 ka.

Deep Pacific and global mean δ 13 C
Previous studies have assumed deep Pacific δ 13 C can be used as a proxy for global mean δ 13 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 δ 13 C stacks show similar patterns of change (Fig. 6) and fall along a tight regression line δ 13 C global = 1.04 ± 0.06 ‰ × δ 13 C DP + 0.19 ± 0.01 ‰ (1) The two time series are highly correlated (r 2 = 0.99), which is not surprising because the large volume of the deep Pacific exerts a strong influence on the global mean δ 13 C stack. When the stacks are pre-whitened to account for autocorrelation (Table 2), their correlation is weaker (r 2 = 0.46) but statistically significant (p = 0.05).
Alternatively, a carbon cycle box model simulated a strong correlation between deep Pacific δ 13 C and CO 2 across several glacial cycles (r 2 = 0.96; Köhler et al., 2010). The correlation between CO 2 and our deep Pacific δ 13 C stack is statistically significant after pre-whitening (r 2 = 0.57, p = 0.02), but global mean δ 13 C and CO 2 are not (r 2 = 0.28, p = 0.16). Our compilation of Pacific records is likely insufficient to determine whether deep Pacific δ 13 C correlates better with global mean δ 13 C or atmospheric CO 2 . This issue could be better resolved using a δ 13 C compilation spanning multiple glacial cycles and including more deep Pacific sites.
6.3 Vertical δ 13 C gradient and atmospheric CO 2 The vertical δ 13 C gradient ( δ 13 C I−D ) in our compilation resembles the inverse of CO 2 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 2 using the gradient between the deep Pacific and half the INA δ 13 C stack ( δ 13 C (INA/2)−DP ), Lisiecki (2010). Both vertical δ 13 C gradients ( δ 13 C I−D and δ 13 C (INA/2)−DP ) decrease from 18-11 ka over the same time interval that atmospheric CO 2 increases. In contrast, the global mean δ 13 C stack increases at a relatively steady pace from 19-6 ka. Thus, the δ 13 C gradients record a distinctly different signal than global mean δ 13 C.
The δ 13 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 2 rise (Marcott et al., 2014;Köhler et al., 2017) and a decrease of 0.3 ‰ in the δ 13 C of atmospheric CO 2 (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 δ 13 C gradients at 12 ka approximately www.clim-past.net/14/1229/2018/ coincides with rapidly increasing atmospheric CO 2 from 13-11.5 ka and a decrease of 0.1 ‰ in the δ 13 C of atmospheric CO 2 .
From 11 to 6 ka, atmospheric CO 2 remains nearly constant with a small (approximately 10 ppm) decrease from 11-8 ka. The vertical δ 13 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 2 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., 2005Köhler et al., , 2010. Although CO 2 correlates strongly with both δ 13 C I−D (r 2 = −0.96) and δ 13 C (INA/2)−DP (r 2 = −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 2 significantly correlates with both δ 13 C I−D (r 2 = −0.51; p = 0.03) and δ 13 C (INA/2)−DP (r 2 = −0.69; p < 0.01) (Fig. S1, Table 2). Better correlation with δ 13 C (INA/2)−DP could be because of better age control and higher resolution δ 13 C records in the INA region than the other intermediate regions. Determining whether the δ 13 C (INA/2)−DP gradient or the global vertical δ 13 C gradient correlates better with atmospheric CO 2 will require data with better spatial coverage and/or a longer time span.
Because our comparison of the vertical δ 13 C gradient and CO 2 could be affected by lags within the carbon cycle or age model uncertainty, we additionally investigate whether the correlations between CO 2 and the vertical δ 13 C gradient would be improved by age model shifts ( Table 2). The correlation between CO 2 is maximized when δ 13 C I−D or δ 13 C (INA/2)−DP lags CO 2 by 400 years (Table 2), which is within the age uncertainty of the sediment core age models. Thus, changes in atmospheric CO 2 and vertical δ 13 C gradients appear synchronous to within age model uncertainty.
Processes that potentially explain atmospheric CO 2 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 2 change. However, because both productivity and circulation change would affect the vertical δ 13 C gradient while changing atmospheric CO 2 , 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 δ 13 C gradients as a proxy for glacial-interglacial CO 2 change on both orbital and millennial timescales (Lisiecki, 2010).

Conclusions
We present regional δ 13 C stacks and volume-weighted intermediate, deep, and global mean δ 13 C stacks from a compilation of 127 benthic C. wuellerstorfi δ 13 C records, which span 20 to 6 kyr with a mean age resolution better than 2 kyr. Age models are based on δ 18 O alignments to regional stacks with radiocarbon dating and age model uncertainties of approximately 1-2 kyr. Our compilation shows spatial patterns in benthic δ 13 C that are similar to higher resolution reconstructions of the Holocene and Last Glacial Maximum. The volume-weighted mean δ 13 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 δ 13 C time series also allows us to evaluate the timing of change in the mean and vertical gradient of δ 13 C and compare them with other carbon cycle changes. The volume-weighted global δ 13 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, δ 13 C changes from 20 to 6 ka suggest that a deep Pacific δ 13 C stack approximates global mean δ 13 C with an offset of 0.19 ‰. Vertical δ 13 C gradients between intermediate and deep water ( δ 13 C I−D ) and between the intermediate North Atlantic and deep Pacific ( δ 13 C (INA/2)−DP ) are interpreted as proxies for change in deep ocean carbon storage. Millennial-scale features in δ 13 C I−D and δ 13 C (INA/2)−DP are significantly correlated with atmospheric CO 2 changes from 20-6 ka.
Based on these analyses, we conclude that the fourdimensional compilation of globally distributed δ 13 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.  Table A1. Supplemental table of the name, location, region, and reference for each record in this δ 13 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 δ 13 C record are provided in Table S1 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 2 record in 100-year increments relative to the δ 13 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. Author contributions. LEL contributed to the design of the experiment, interpretation of results, and manuscript revision. CDP contributed data collection, performed the analysis, interpreted results, and wrote the manuscript.
Competing interests. The authors declare that they have no conflict of interest.

Special issue statement.
This article is part of the special issue "Paleoclimate data synthesis and analysis of associated uncertainty (BG/CP/ESSD inter-journal SI)". It is not associated with a conference.