A synthesis of marine sediment core δ 13 C data over the last 150 000 years

. The isotopic composition of carbon, δ 13 C, in seawater is used in reconstructions of ocean circulation, marine productivity, air-sea gas exchange, and biosphere carbon storage. Here, a synthesis of δ 13 C measurements taken from foraminifera in marine sediment cores over the last 150 000 years is presented. The dataset comprises previously published and unpublished data from benthic and planktonic records throughout the global ocean. Data are placed on a common δ 18 O age scale suitable for examining orbital timescale variability but not millennial events, which are removed by a 10 ka ﬁlter. Error estimates account for the resolution and scatter of the original data, and uncertainty in the relationship between δ 13 C of calcite and of dissolved inorganic carbon (DIC) in seawater. This will assist comparison with δ 13 C of DIC output from models, which can be further improved using model outputs such as temperature, DIC


Abstract.
The isotopic composition of carbon, δ 13 C, in seawater is used in reconstructions of ocean circulation, marine productivity, air-sea gas exchange, and biosphere carbon storage. Here, a synthesis of δ 13 C measurements taken from foraminifera in marine sediment cores over the last 150 000 years is presented. The dataset comprises previously published and unpublished data from benthic and planktonic records throughout the global ocean. Data are placed on a common δ 18 O age scale suitable for examining orbital timescale variability but not millennial events, which are removed by a 10 ka filter. Error estimates account for the resolution and scatter of the original data, and uncertainty in the relationship between δ 13 C of calcite and of dissolved inorganic carbon (DIC) in seawater. This will assist comparison with δ 13 C of DIC output from models, which can be further improved using model outputs such as temperature, DIC concentration, and alkalinity to improve estimates of fractionation during calcite formation.
High global deep ocean δ 13 C, indicating isotopically heavy carbon, is obtained during Marine Isotope Stages (MIS) 1, 3, 5a, c and e, and low δ 13 C during MIS 2, 4 and 6, which are temperature minima, with larger amplitude variability in the Atlantic Ocean than the Pacific Ocean. This is likely to result from changes in biosphere carbon storage, modulated by changes in ocean circulation, productivity, and air-sea gas exchange. The North Atlantic vertical δ 13 C gradient is greater during temperature minima than temperature maxima, attributed to changes in the spatial extent of Atlantic source waters. There are insufficient data from shallower than 2500 m to obtain a coherent pattern in other Correspondence to: K. I. C. Oliver (k.oliver@noc.soton.ac.uk) ocean basins. The data synthesis indicates that basin-scale δ 13 C during the last interglacial (MIS 5e) is not clearly distinguishable from the Holocene (MIS 1) or from MIS 5a and 5c, despite significant differences in ice volume and atmospheric CO 2 concentration during these intervals. Similarly, MIS 6 is only distinguishable from MIS 2 or 4 due to globally lower δ 13 C values both in benthic and planktonic data. This result is obtained despite individual records showing differences between these intervals, indicating that care must be used in interpreting large scale signals from a small number of records.

Introduction
The isotopic composition, δ 13 C, of inorganic carbon in seawater is a diagnostic of ocean circulation and the marine and terrestrial carbon cycle. The potential for δ 13 C in calcium carbonate shells, formed by foraminifera and preserved in marine sediments, to record past changes in climate has been recognised since the 1970s (Shackleton, 1977;Duplessy et al., 1981a). Greater differences in δ 13 C between planktonic and benthic foraminifera, during glacial periods, were interpreted as indicating greater storage of isotopically light organic carbon in the deep ocean and linked to atmospheric pCO 2 (Broecker, 1982a;Shackleton et al., 1983;Shackleton et al., 1992). Lower δ 13 C values recorded in the Pacific Ocean, during the last glacial maximum (LGM), were attributed to a change in the ocean δ 13 C reservoir due to the release of organic carbon from the terrestrial biosphere or marine shelf sediments (Broecker, 1982b;Duplessy et al., 1988b). Carbon isotopes have also been used to reconstruct past water masses, notably low δ 13 C Antarctic Light microenvironment formed around infaunal foraminifera; may also form around epifaunal foraminifera under high productivity conditions (phyto-detritus effect) Secondary calcification on planktonic species Fig. 1. Schematic of processes influencing δ 13 C recorded in calcium carbonate from foraminiferal shells. Processes in large font affect δ 13 C of dissolved inorganic carbon (DIC) in open seawater; processes in small font affect the difference between recorded δ 13 C and seawater δ 13 C. POC ≡ particulate organic carbon. Heavy (light) DIC refers to DIC rich (poor) in 13 C.
source waters and high δ 13 C Atlantic source waters (Sarnthein et al., 1994), or to diagnose changes in air-sea gas exchange (Marchitto and Broecker, 2006). However, the interpretation of the δ 13 C is complicated by the dependence of fractionation during calcification on properties of the water (temperature; [CO 2 ]; [CO 2− 3 ]) and the foraminifer (species; shell size; diet), and on the formation of microenvironments around benthic foraminifera (Sect. 3). The variety of mechanisms influencing the δ 13 C record is summarised in Fig. 1.
The reconstruction of past ocean states depends on the collation of δ 13 C observations from throughout the ocean, which can provide a fairer test to hypotheses than individual datasets. The demand for such syntheses is increased by the development of Earth system models that are able to simulate paleoclimate proxies including δ 13 C (Ridgwell et al., 2007;Brovkin et al., 2007). Previous data syntheses have consisted of time-slices for selected regions such as the Atlantic Ocean (Sarnthein et al., 1994;Bickert and Mackensen, 2003;Curry and Oppo, 2005;Lynch-Stieglitz et al., 2007), typically focusing on the difference between the late Holocene and the LGM. Data from the benthic species Cibicidoides wuellerstorfi, considered the most reliable indicator of seawater δ 13 C, were selected for those time-slices, reducing the error in the data at the expense of reduced data coverage. In this study, we focus not on specific time-slices but on producing a synthesis of time-slices for the last 150 ka. We include data from a range of benthic and planktonic species, and do not exclude data from high productivity regions. The inhomogeneity in this data is addressed by attaching an error estimate for each observation, using the entire dataset to estimate additional errors associated with less reliable species and unfavourable core locations.
The goals for this study are threefold: (1) to provide a common δ 18 O-derived age-scale for a synthesis of δ 13 C data; (2) to provide δ 13 C error estimates for the synthesis, in order to facilitate direct model-data comparison; (3) to provide a global overview of the data, and an account of the large scale processes invoked to explain changes in ocean δ 13 C. In Sect. 2, we introduce the data and the age-scale. In Sect. 3, we describe the methods used to determine uniformly spaced time-series and uncertainty intervals. In Sect. 4, we examine δ 13 C time-series, grouped by region, and time-slices, as well as planktonic-benthic differences. We summarise our findings in Sect. 5, and discuss their application for Earth System models and biosphere reconstructions. Table 1 summarises the data used in this compilation. Data consist mostly of records submitted to the PANGAEA publishing network for geoscientific and environmental data (http://pangaea.de), the National Geophysical Data Centre (NGDC; http://www.ngdc.noaa.gov), or the Delphi Project (http://www.esc.cam.ac.uk/research/ research-groups/delphi), with additional records obtained through personal correspondence. δ 13 C records were accepted into the pre-processed compilation provided that it was possible to obtain a reliable (though sometimes low resolution) age model within the last 150 000 years using δ 18 O stratigraphy. No data quality constraint was applied to the pre-processed compilation, but several factors (detailed in Sect. 3) can lead to a large error estimate. Where the cumulative errors exceeded 1‰, data were rejected as having too large an error to be useful. Entire records rejected by this process are flagged in Table 1. Figure 2 shows the core locations for the records used in the data synthesis, as well as for records for which we were unable to provide a δ 18 O-derived age model. Among the records that were used, there is good data coverage in deep waters (>2500 m) in the Atlantic Ocean (100 cores with benthic data). Atlantic data coverage is also reasonable at shallower depths (45 cores) and in surface waters (71 cores with planktonic data). The principal gaps in benthic data coverage are: (1) <2500 m in the Pacific ocean (10 cores); (2) the Indian ocean (12 cores); and (4) south of 50 • S (2 cores). There is a large variation in the temporal resolution of the pre-processed data, both within and across cores. Of the cores providing benthic data, all but 28 provide at least one observation for the LGM (19-23 ka), compared with 64 cores providing no Holocene (<10 ka) data, and 124 cores providing no data from Marine Isotope Stage (MIS) 5a or earlier (>71 ka). This is summarised in Table 1. The raw data and the δ 18 O age model are included as Supplementary Materials. The raw data has undergone basic filtering; for example, Table 1. Cores used in the data synthesis. Flags: R = mean data resolution poorer than 6 ka; r = mean data resolution poorer than 3 ka; H = at least one period of >20 ka with no obseravtions; h = at least one period of >10 ka with no observations; p = upwelling region or phytodetritus effect suspected; X = no processed data from record. Abbreviations are as follows. In header rows: Plank = planktonic; # of obs. = number of observations in core; Sed. rate = mean sedimentation rate (cm/ka) over period of data; all latitudes and longitudes in decimal degrees. In Ocean column: Arct. = Arctic Ocean; Arct N = Nordic seas; C = Caribbean Sea; S = Southern Ocean (defined as >40 • S); SC = South China Sea. In Flag columns; R = mean data resolution poorer than 6 ka; r = mean data resolution poorer than 3 ka; H = at least one period of >20 ka with no obseravtions; h = at least one period of >10 ka with no observations; p = upwelling region or phytodetritus effect suspected; X = no processed data from record. Relevant references are as follows:  where records consisted principally of data from a single species, but with additional observations from other species, these additional observations were removed. We caution that this is not a primary source of raw data, which should instead be obtained from PANGAEA, NGDC, or Delphi. The fully processed data, after applying the methods described in Sect. 3, are included as Supplementary Materials.

Age modelling
For this study it was necessary to construct consistent age models for all cores. Since the accuracy required of these age models is rather low (a good representation of precessional and longer timescale variability is sufficient), we maximised the number of records that could be included by constructing age models based solely on oxygen isotope data. A standardised framework was developed that allows datasets to be placed on a common chronological framework for the last 150 000 years. For this we determined the timing of δ 18 O maxima that can be recognized in both benthic and planktonic δ 18 O records, using the chronology of Shackleton et al. (2000), which can be readily recognized in higher resolution records. Approximately half of the benthic oxygen isotope signal results from changes in sea-level that affect benthic δ 18 O synchronously to within the mixing timescale of the ocean. Therefore, δ 18 O maxima approximate to sea level minima (Shackleton et al., 2000). The ages of δ 18 O maxima -18, 62, 87, 108 and 137 ka -are compatible with the ages of analogous events at Dome Fuji on a timescale based on the O 2 /N 2 ratio of gas trapped in the ice-core (Kawamura et al., 2007), which has 2σ absolute age uncertainties of approximately 1000-2500 years. Dates between δ 18 O maxima were estimated by applying uniform sedimentation rates between δ 18 O maxima. The resulting δ 18 O timeseries, filtered by the same method as δ 13 C data (described in Sect. 3.1.1) are plotted in Fig. 3, as anomalies relative to the LGM. δ 18 O maxima were initially treated as synchronous throughout the ocean, despite known asynchroneities in climate events expected to distort the δ 18 O signal. Asynchroneities of 1.5-3 ka have been identified between the Northern and Southern ice sheets (Blunier and Brooks, 2001). Using 14 C dating, Skinner and Shackleton (2005) and Labeyrie et al. (2005) obtained Atlantic-Pacific leads of 3.9 and 1.5 ka, respectively, for the most rapid decrease in δ 18 O during deglaciation, and Labeyrie et al. (2005) obtained large asynchroneities between intermediate and deep waters within the Atlantic basin. Due to the sparsity of 14 C records, and their inapplicability to the 62, 87, 108, and 137 ka δ 18 O maxima, systematic regional changes to the age models were not made. Moreover, both the Skinner and Shackleton (2005) and Labeyrie et al. (2005) studies, δ 18 O maxima from all records were within 2000 years of 18 ka, which we show to be small compared with other likely errors in the age models.
We now provide an estimate of the age error associated with our approach and discuss proposed improvements. Age errors in individual records could be reduced using additional stratigraphic constraints or age control points from 14 C data or tephra layers. However, such data were not available for most records, and to introduce heterogeneity to the synthesis by applying a mixture of methods would likely increase the relative age errors between cores. Instead, we estimate the   error associated with our age modelling approach by comparing our age estimates with previously published age models for the 116 cores for which these were available with the raw data (see Table 1 for references). Figure 4 shows the distribution of offsets of our age estimates from previous age estimates. For each of the four ocean basins considered, the age models agree to within 6 ka for approximately 95% of observations. Since there are more observations from high resolution records, this result is biased towards such records, whereas greater errors might be expected in low resolution records. Giving each record equal weight by normalising by the record length, 93% of observations agree to within 6 ka. The distribution is approximately gaussian, with no clear bias towards underestimating or overestimating age in any basin. This is also true if the previous age estimates used are restricted to those based on 14 C measurements (plotted in bold in Fig. 2). In Sect. 4.4 we show that applying regionally uniform adjustments to our age model in order to minimise the age error in records with 14 C derived age models has negli-gible effect on our results, and therefore we persist with the initial approach of treating δ 18 O maxima as synchronous.
In ascribing age uncertainties, we note that disagreement between two age models may indicate error in either age model, and may be associated with the error in the absolute age scale rather than in the relative ages of records, suggesting that it would be pessimistic to equate the difference between the two age estimates with the error in our age estimates. On the other hand, errors common to our age estimates and previously published estimates are possible where 14 C data were not used in constructing either age model (for example, asynchroneities between δ 18 O signals are likely to affect orbitally tuned age models similarly to our age model). We ascribe a 2σ error of 6 ka throughout most of the dataset, slightly greater than the error of 5 ka attributed to age models obtained by orbitally tuning of several proxies by Martinson et al. (1987). Errors are much smaller between the LGM and the Holocene, and typically greatest during Marine Isotope Stage (MIS) 3 (29-57 ka), which we attribute  (Table 1) age estimates as a function of age. Colours indicate ocean basin; bold filled circles indicate that 14 C data was used in the previously published age model (note that 14 C can only be used to improve age estimates for the last ∼50 ka). The bold black line indicates offsets of ±6 ka. (b) Distribution of age estimate offsets. (c) Distribution of age estimate offsets, normalised by the number of observations within each record (i.e. an observation from a record of length n adds 1/n, as opposed to 1, to the appropriate bin).
to the large interval (between 18 and 62 ka) over which a uniform sediment accumulation rate must be assumed, and MIS 6 (>137 ka). The ascribed age error is increased to 8 ka during these MIS 3 and 6. We attribute the increased errors during MIS 3 to the long interval between δ 18 O maxima (44 ka), which would be expected to increase errors due to non-uniform sedimentation rate. We therefore emphasise that, due to the age modelling approach employed, the synthesis is not suitable for examining the phasing of climate events between ocean basins, nor millennial timescale variations in δ 13 C.

Procedure for temporal gridding and δ 1C uncertainty estimation
This section describes the treatment of data inhomogeneity within and between records in obtaining a temporally gridded data synthesis, with uncertainties ascribed to individual δ 13 C values as a function of species, core location, data resolution and scatter. In accounting for data inhomogeneities, it is useful to distinguish between systematic errors associated with each record (inter-record error) and the temporally varying error component (intra-record error). For example, the error due to fractionation during calcification has a large inter-record component because this is strongly dependent on species, and most records consist of samples from a single species or genus. Spatial variability in the species offset is also an inter-record error, but in situ variability of this offset in response to climate change is an intra-record error. Similarly, any systematic difference in δ 13 C resulting from inter-laboratory calibration errors is an inter-record error. Since we expect uncertainties in the systematic inter-record error to be significant, δ 13 C and its uncertainties are estimated by two stages. First, we present intra-record variability in δ 13 C relative to the LGM, which we label LGM δ 13 C (Sect. 3.1). Second, we obtain estimates for absolute δ 13 C values (Sect. 3.2).
LGM δ 13 C can be used to understand spatial patterns of δ 13 C variability, and has smaller uncertainty than absolute δ 13 C. Therefore, we recommend that LGM δ 13 C is chosen where possible when using this data synthesis (for example for model-data comparison), and it is also the variable presented in most figures within this paper.

Obtaining δ 13 C as anomalies relative to the LGM
A smoothing spline, approximating a 10 ka low pass filter (Sect. 3.1.1), was applied to the data, removing high frequency variability such as Dansgaard-Oeschger oscillations which are beyond the precision of the age model. All LGM δ 13 C values are calculated as an anomaly relative to the LGM. This is defined at 21 ka, providing a good approximation to the LGM definition of 19-23 ka (e.g. Mix et al., 2001); averaging pre-processed data within the 19-23 ka interval would not be preferable due to uncertainties in the age model. The LGM, rather than the late Holocene, is chosen as the reference period because almost all cores provide LGM data, whereas approximately half of all cores provide no data for < 4 ka.
Three steps were taken to obtain LGM δ 13 C estimates and uncertainties, described in detail below. First, each timeseries was smoothed and gridded, and the error associated with this process estimated (Sect. 3.1.1). Second, an additional error was estimated for cores thought to be subject to changes in the phyto-detritus effect (Sect. 3.1.2). Third, an error for the representativeness of each species of the ambient water was estimated (Sect. 3.1.3). Although the errors associated with each of these steps are likely to be correlated, there were insufficient data to reliably estimate these correlations, so the total error was estimated to be the quadratic sum of the error components. 2σ error estimates are quoted; where they exceed 1‰ or where the data resolution is poorer than 6 ka, no final LGM δ 13 C estimate is given. Full error estimates are provided as part of the uniformly spaced output dataset (Supplementary Materials). Data are output on a 2-ka grid as an aid to visualisation, but note that the age error estimates are of 6-8 ka.

Smoothing splines
The first step in obtaining LGM δ 13 C estimates was independent of species or region, and a function only of the input time-series data. These were smoothed and placed on a uniform 2 ka grid using the spline smoothing method described by Silverman (1985). For a time-series of n observations, the set of n − 1 cubic equations was obtained which minimised where λ is the smoothing parameter, y is the output δ 13 C estimate andÿ is its second derivative with respect to time, t 1 and t n are dates for the top and bottom, respectively, of each estimate offsets, normalised by the number of observations within each record (i.e. an observation from of length n adds 1/n, as opposed to 1, to the appropriate bin).  record (note that t is defined positive towards the past), and y is given by where Y i is the δ 13 C observation and y thr = 0.5‰ was chosen to prevent outliers from having a dominant effect on the estimate (Enting, 1987). The smoothing parameter is given by where T 0.5 was chosen to be 10 ka. For a uniformly spaced input time-series, this method is equivalent to a kernel filter admitting 50% of variability with a period of T 0.5 , steeply increasing towards 100% at longer timescales and steeply decreasing towards 0% at shorter timescales (Enting, 1987). For time-series with non-uniformly spaced data, the period at which 50% of the variability is admitted is a weak function (∼f − 1 4 ) of the local kernel densityf . Kernel data density estimates, which approximate to the inverse of the temporally local data resolution t, are provided as part of the uniformly spaced output dataset (Supplementary Materials).
A first estimate of the error in the smoothed time-series is given by (Silverman, 1985) 2σ Typical values of 2σ sil are of the order of 0.1‰, rarely increasing above 0.25‰. However, this method is intended for time-series with large n, and unrealistically small error estimates are produced for time-series with a temporal resolution of the same order as T 0.5 . Enting (1987) noted that an accurate error estimate depends on a knowledge of the spectral characteristics of both the signal and the noise. For our purposes, all high frequency variability is noise, and it is likely that this is aliased in low resolution time-series. In order to address this, we assume that the root-mean-square of the residuals,ȳ , from high resolution time-series provides an estimate of the noise present in all time-series. Where smaller values ofȳ are obtained from low resolution records, this indicates that variability in the high frequency (T < T 0.5 ) band is being aliased as a signal at lower frequencies.ȳ is plotted as a function of time-series mean resolution in Fig. 5, along with a line of best fit, y h , with the form of a tanh function. We obtain lower values of y h in low resolution cores. The resulting error estimate due to aliasing, is also plotted in Fig. 5. Note that there is a large range in y h even within high frequency time-series, likely to be largely caused by data inhomogeneities about which information is often unavailable (e.g. shell sizes used; shells per sample; averaging of duplicate samples). Therefore aliasing error estimates for low resolution time-series are very uncertain. The total error associated with the smoothing spline process is the quadratic sum of 2σ sil and 2σ ali .

Error due to the phyto-detritus effect
A further source of error is the possible change in time of the depositional environment in which benthic foraminifera exist. Infaunal species such as those in the Uvigerina genus inhabit a microenvironment that is often depleted in δ 13 C due to the respiration of organic matter in the sediment, whereas epifaunal species are not typically subject to this effect. There is some evidence that the offset in infaunal species varies over a glacial cycle (Hoogakker et al., 2010), but there are insufficient data to determine a correction that varies over time. In this study, we treat this effect as a constant component of the species offset plus a random error (Sect. 3.1.3); therefore there is no phyto-detritus correction to LGM δ 13 C in infaunal species. Although shells from the epifaunal Cibicidoides taxon are thought usually to record the δ 13 C of ambient water CO 2 accurately, several studies have identified a bias toward low values in high productivity regions (Mackensen et al., 1993;Sarnthein et al., 1994), termed the phyto-detritus effect. Bickert and Wefer (1999) compared cores from upwelling regions in the Atlantic Ocean to nearby cores outside upwelling regions and found that δ 13 C in upwelling regions was more depressed, relative to nearby cores, during the LGM than the late Holocene. Bickert and Mackensen (2003) applied a correction of +0.4‰ to LGM δ 13 C in several upwelling regions in the Atlantic Ocean (parts of the eastern boundary in both hemispheres; equator; subtropical and subantarctic fronts), applying no analogous correction to Holocene values. Such a time-dependent correction would influence both δ 13 C and LGM δ 13 C estimates. However, the temporal and spatial structure of such a correction is poorly understood. Instead, we add a 2σ phy error of 0.4‰ to the records upwelling regions where phytodetritus errors are more likely (flagged in Table 1). In Sect. 4.4 we investigate whether a bias is introduced to the data synthesis by not applying a phytodetritus correction.
It is not well known to what extent the phyto-detritus effect in epifaunal species is caused by changes in foraminifer growth rate (McConnaughey et al., 1997), or by the build-up of respiring organic matter over epifaunal species leading to a decrease in δ 13 C, similar to that observed in infaunal species. No correlation has been observed between sedimentation rate and amplitude of the phyto-detritus effect, although Mackensen et al. (2001) suggested that this is because the key control is seasonal peak deposition rate, rather than annual mean deposition rate. It is therefore uncertain whether additional errors are introduced to records from infaunal species in high productivity regions. As noted above, we apply no correction, but add the same uncertainty as would be added to epifaunal species.

Disequilibria and the calculation of species errors
An important error in LGM δ 13 C estimates is variability in disequilibria of δ 13 C in calcium carbonate in the shells of foraminifera, relative to the ambient water. Along with the phyto-detritus effect, δ 13 C disequilibrium leads to the application of different "species offsets" to δ 13 C data collected from different species. We note that LGM δ 13 C is independent of any uniform species offset, and that errors due to variability between individual shells, or small groups of shells, of a given species have already been implicitly accounted for in the smoothing spline calculation. Of interest here are the causes of variability in the disequilibrium of calcium carbonate through time.
In planktonic formanifiera, recorded δ 13 C is a strong function of shell size (e.g. Berger et al., 1978;Spero and Lea, 1993;Elderfield et al., 2002), and the difference in δ 13 C between the largest and smallest size fractions can be greater than glacial-interglacial differences of the order of 0.5‰ (Oppo and Fairbanks, 1989). Therefore, δ 13 C measurements from both planktonic and benthic species are usually made on shells selected for size fraction. Nevertheless, any unrecorded changes in shell size throughout a record may lead to significant error in LGM δ 13 C. A second source of error is temperature-dependent fractionation, which affects planktonic species due to variability in surface temperature on seasonal to glacial timescales. The temperature-dependence of fractionation is a function of species. Culture experiment estimates of ∂(δ 13 C)/∂(T ) (where T is temperature) have been made of −0.11‰ • C −1 for Globigerina bulloides (Bemis et al., 2000), −0.08‰ • C −1 for Limacina inflata , 0 to +0.05‰ • C −1 for Orbulina universa (Bemis et al., 2000). No relationship with temperature was Table 2. Species-specific error estimates and uniform adjustments applied to the dataset. Correction and correction error estimates affect the calculation of δ 13 C but not LGM δ 13 C; where no correction is given, no absolute δ 13 C estimate is made.

Species
LGM δ 13 C error, ‰ Correction, ‰ Correction error, ‰  , whereas Kohfeld et al. (2000) assumed a relationship of Culture experiments on planktonic foraminifera also reveal a strong dependence on the concentration of the carbonate ion, with a gradient of approximately −0.012‰/(µmol kg −1 ) for G. bulloides under constant DIC (Spero et al., 1997). Glacial carbonate concentration is dependent on the poorly constrained alkalinity inventory of the ocean, but a change of 30 to 80 µmol kg −1 is not implausible (Kohfeld et al., 2000), indicating an effect of the same order as glacial δ 13 C cycles. The effect of carbonate ion on benthic δ 13 C is not known. Increases in light availability have been found to increase δ 13 C in Globigerinoides sacculifer, leading Spero and Lea (1993) to propose that the correlation between shell size and δ 13 C is principally caused by the dependence of both growth rate and δ 13 C on ambient light levels. Further effects arise due to changes in the isotopic composition of carbon in the diet. Culture experiments indicate that δ 13 C in G. bulloides shell carbonate varies as δ 13 C in organic matter in the ratio 0.084:1 (Spero and Lea, 1996). Finally, the planktonic δ 13 C signal can be distorted from a surface signal by secondary calcification in deeper, typically isotopically lighter, waters, as observed by Duplessy et al. (1981b) on G. sacculifer transferred to laboratory tanks.
With the exception of sea-surface temperature reconstructions for the LGM (e.g. MARGO Project Members, 2009), variability in water properties likely to cause systematic errors in LGM δ 13 C are poorly constrained over the last 150 ka, so no time-varying correction is applied to LGM δ 13 C. Therefore, we caution that changes in these properties over a glacial cycle, rather than changes in seawater δ 13 C, may contribute significantly to variability in the planktonic LGM δ 13 C record. However, Earth system models that output δ 13 C can also produce temperature, alkalinity and pH fields, so that this is not an obstacle to model-data comparison; the validity of such comparisons are instead limited by the general applicability of the above empirical relationships. Moreover, we note that these caveats are a much lesser concern when interpreting the benthic δ 13 C record.
A random error estimate may be made by comparing records from different species but from the same core. Among benthic species, the epibenthic taxon Cibicidoides wuellerstorfi is widely preferred for recording δ 13 C, and to this species we ascribe an error of 0.15‰. This value is obtained from comparison of water column data (Kroopnick, 1985) with data from core samples (c.f. Beveridge et al., 1995). Error estimates for other species in the Cibicidoides genus were obtained using cores containing both C. wuellerstorfi and another Cibicidoides species. The error is the rootmean-squared (rms) difference between the two species after the application of an optimal species offset. This process was repeated for the infaunal genus Uvigerina. For other species, there was insufficient overlap in data with C. wuellerstorfi to apply this method. For these species, errors of 0.4‰ were ascribed. A larger error of 0.6‰ was applied to data derived from Hoeglandia elegans or from mixed benthics. Error estimates are presented in Table 2.
For planktonic foraminifera, the choice of reference species is less clear. We select G. ruber, which yields data that are consistent on a regional scale, and has a restricted depth range , so that changes in recorded δ 13 C are more likely to reflect changes in surface conditions. We ascribe an error of 0.4‰ to this species, and calculated estimates for G. bulloides, Globigerinoides sacculifer, and N. pachyderma by the same method used for benthic species. However, since these relative errors are less than 0.4‰, we ascribe an error of 0.4‰ to each of these species. These are presented in Table 2. We ascribe an error of 0.6‰ to other planktonic species.

Obtaining absolute δ 13 C values
Absolute δ 13 C values within each time-series are the sum of LGM δ 13 C estimates and the LGM δ 13 C value. The LGM δ 13 C estimate is the value output from spline smoothing plus a species-specific uniform correction. For the benthic genus Cibicidoides and the planktonic species G. ruber, we ascribe a correction of 0 ± 0.2‰. The correction and correction error for other species are obtained by calculating the mean and standard deviation of the offsets used to optimise the leastsquares-fit between different species in the same core, as described in Sect. 3.1.3. This precludes any spatial variability in the species offset, because despite evidence that these offsets change over a glacial cycle (Hoogakker et al., 2010), there are insufficient data to quantify this variability globally. We obtain a correction of +0.85 ± 0.48‰ for Uvigerina species, similar to the canonical value of +0.9‰ (Shackleton and Hall, 1984). Other estimates are given in Table 2; where no correction is quoted, no absolute δ 13 C estimate is made.

δ 13 C data presentation and interpretation
For the purpose of visualisation, data were sorted into the principal regions: North Atlantic, South Atlantic, Indian Ocean and Pacific Ocean (Table 1), with subcategories for regions such as the Arctic Ocean, the Nordic Seas, and the South China Sea, and the Southern Ocean sector of each ocean, for which a highly inclusive definition of south of 40 • S is used. By compiling all available δ 13 C data in each region, presented as an anomaly relative to the LGM, and plotting on a uniform timescale, we can look for large scale changes in δ 13 C that might be obscured or biased by consideration of a small number of cores. Nevertheless, there are sampling biases towards coastal areas in each ocean, towards the eastern Atlantic and towards the Arabian Sea in the Indian Ocean (Fig. 2), as well as towards the 2500-3500 m depth range. Data coverage is insufficient to interpolate and extrapolate in three dimensions, but estimates might be refined with new data or by dynamical smoothing using an Earth system model.   (Kawamura et al., 2007); note that the age model used in this study is based on the Kawamura et al. (2007) timescale.

Benthic time-series
In Fig. 6 we present LGM δ 13 C time-series grouped by region, excluding all data from the Arctic Ocean or from marginal seas (except the South China Sea), and by depth.
Only where the temporal resolution of the pre-processed data was sufficiently high ( t < 6 ka) and the estimated error suf-ficiently low (<0.8‰) are data plotted. The mean and 2× standard error of data from each region were calculated and are plotted in Fig. 7. Sampling biases towards ocean boundaries, and towards regions with high sedimentation rates, are not removed in the averaging calculation. Therefore Fig. 7 is not an unbiased estimate of regional averages, but rather an estimate of LGM δ 13 C where data exist within each region. In waters deeper than 2500 m, a consistent pattern is observed within the North Atlantic, South Atlantic, and Pacific Oceans. These basins exhibit maxima in δ 13 C during temperature maxima, with small differences between the Eemian interglacial (MIS 5e), the MIS 5c and MIS 5a interstadials, and the Holocene (MIS 1). MIS 4 and MIS 6 are marked by values of δ 13 C similar to or lower than LGM values, consistent with a longer timescale trend towards higher δ 13 C (Hoogakker et al., 2006). The principal difference between the basins is in the amplitude of δ 13 C variations. The North Atlantic exhibits slightly larger amplitude changes in δ 13 C than the South Atlantic. The amplitude of the glacial cycle in the Pacific basin is approximately half that in the Atlantic basin. Furthermore, there is much less variability observed in the Pacific at around precessional timescales; MIS 3 is a much weaker peak, whereas troughs during MIS 5b and MIS 5d are not consistently recorded. Typically lower sedimentation rates in the Pacific, leading to poorer data resolution, may contribute to the latter result. Another bias may arise from a weak positive correlation between bathymetric depth and the amplitude of variability. However, the mean depth of >2500 m cores in the Pacific (3400 m) is only slightly shallower than that in the Atlantic (3600 m), so this can explain less than 0.05‰ of the difference between the two basins. The lack of long records from deeper than 2500 m in the Indian Ocean prevent us from concluding whether there is similar variability in δ 13 C over a full glacial cycle, although the amplitude of change between the LGM and the Holocene is greater than typical Pacific but smaller than typical Atlantic amplitudes.
Much of the variability in deep LGM δ 13 C is likely to arise from changes in the global ocean δ 13 C reservoir, influenced by the storage of isotopically light carbon in the biosphere. This is supported by the qualitative similarity between the Pacific and the Atlantic LGM δ 13 C, unless there are changes of the opposite sign within intermediate or pycnocline waters. There is some evidence, from Pacific cores in the 1500-2500 m depth range, for such changes over the time interval since the LGM (c.f. Duplessy et al., 1988b). However, records from a similar depth in the Indian Ocean, as well as shallower Pacific records, present a pattern that is broadly consistent with that observed in the deep Pacific.
Glacial variations in intermediate and pycnocline Atlantic δ 13 C are thought to be strongly influenced by changes in the depth and range of the Atlantic overturning cell, in competition with Antarctic Bottom Water and Antarctic Intermediate Water, although changes in the partitioning of remineralised carbon have also been invoked (Boyle, 1988). Here, we note that Atlantic LGM δ 13 C in the 1500-2500 m depth range is similar to that observed in the Indian Ocean at the same depth, and in the deep Pacific Ocean. An exception is during deglaciations, when there are lower values of Atlantic intermediate water LGM δ 13 C. Holocene LGM δ 13 C in <1500 m deep cores is usually negative, particularly in the West Atlantic, whereas there is no consistent pattern prior to the LGM. ig. 9. Mean and 2× standard error, weighted by the square of the inverse of the error estimate, of planktoni inus benthic ∆ LGM δ 13 C (∆ LGM δ 13 C p−b )in each region. Note that ∆ LGM δ 13 C p−b was calculated for eac cord prior to averaging. Fig. 9. Mean and 2× standard error, weighted by the square of the inverse of the error estimate, of planktonic minus benthic LGM δ 13 C ( LGM δ 13 C p−b )in each region. Note that LGM δ 13 C p−b was calculated for each record prior to averaging.

Planktonic time-series and planktonic/benthic differences
Planktonic LGM δ 13 C time-series, grouped by region, are presented in Fig. 8. There is less agreement between different planktonic records within a region, compared to benthic records, consistent with the large uncertainty attributed to the estimation of seawater LGM δ 13 C from planktonic species (Sect. 3.1.3). G. ruber records from the North Atlantic and Indian Oceans consistently show elevated Holocene LGM δ 13 C values. Temperature maxima also exhibit higher LGM δ 13 C values than temperature minima in most G. bulloides records, although these changes are rarely outside the formal uncertainty limits. Notable exceptions are three cores from the Pacific sector of the Southern Ocean, in which low early Holocene LGM δ 13 C values are recorded, and Indian Ocean cores, which exhibit little coherent variability over 150 000 years, except for a slight upward trend. A decrease in LGM δ 13 C during MIS 4 is present in most Atlantic records but less apparent in Pacific records. There is no clear evidence for minima in planktonic LGM δ 13 C during MIS 5b and MIS 5d, in contrast to benthic records from the deep Atlantic Ocean.
The difference between δ 13 C recorded in planktonic and benthic foraminifera has long been of interest to paleooceanographers because of its association with deep ocean carbon storage (Broecker, 1982a;Shackleton et al., 1983), although interpreting changes in seawater δ 13 C from planktonic data is problematic (Sect. 3.1.3). The difference between planktonic and benthic LGM δ 13 C, LGM δ 13 C p−b , was calculated for each core where both planktonic and benthic δ 13 C records were present. The weighted mean of LGM δ 13 C p−b was calculated for each region and is plotted in Fig. 9 (analogous to Fig. 7). There are few benthic/planktonic pairs from the Indian Ocean, or from shallower than 2500 m in the Pacific Ocean, so that averages are constructed from 3-8 benthic/planktonic pairs in each of these regions (averages are not constructed from fewer than three pairs). Throughout all regions, LGM δ 13 C p−b remains within two standard errors of zero for most of the 150 ka record. The most notable features occur during deglaciations (MIS 1/2 boundary; MIS 5e/6 boundary), where LGM δ 13 C p−b is usually negative outside the North Atlantic. The very low LGM δ 13 C p−b values obtained for MIS 1 in the Pacific are caused by anomalously low planktonic LGM δ 13 C in the Southern Ocean (Fig. 8); exclusion of these cores removes this feature.
Changes in the difference in δ 13 C between surface and deep waters could be caused by changes in the relative strength of export production and meridional overturning in the global ocean. This is because the positive vertical δ 13 C gradient, set up by the downward flux of isotopically light organic matter, is eroded by the exchange of upper ocean and deep waters. The weak δ 13 C gradient found during deglaciations (i.e. the low value of LGM δ 13 C p−b ) may indicate a decrease in export production, more rapid mixing between the deep and surface ocean, or a shift in downwelling towards isotopically heavier parts of the surface ocean. However, several other mechanisms influence LGM δ 13 C p−b , such as regional variability in productivity, circulation, and air-sea gas exchange, and periods of rapid growth or decline of the biosphere. It is also possible that variability in the LGM δ 13 C p−b is dominated by variability in isotope disequilibria in planktonic foraminifera due to changes in temperature, diet, or carbonate ion concentration (Kohfeld et al., 2000).

Time-slices
A set of LGM δ 13 C time-slices, as a function of depth within each ocean basin, for the Holocene (5 ka), MIS 3 (47 ka), MIS 5a (75 ka), and MIS 5e (123 ka) are plotted in Figs. 11-14. These are referenced to an absolute δ 13 C time-slice for the LGM (21 ka), plotted in Fig. 10. We first consider the hypothesis that benthic LGM δ 13 C in each time-slice is a function of ocean region and depth only, and that all deviations from the mean LGM δ 13 C profile can be explained by the quoted error estimates. If the random errors are underestimated in this study, or if LGM δ 13 C varies sufficiently between different water masses at the same depth in the same region to be distinguishable from random error, the apparent scatter in the data would exceed that predicted by the error bars. Since the scatter in the data in Figs. 10-14 are consistent with the quoted errors, we conclude that we have not underestimated the random error in the data. Furthermore, the estimated errors are likely to be too great to use the synthesis to reconstruct local changes in water mass boundaries, for which abundant well dated records from C. wuellerstorfi would be needed. The following discussion emphasises basin-scale changes in δ 13 C that are coherent and statistically significant.
LGM δ 13 C decreases with depth throughout the Atlantic Ocean (Fig. 10), with steep gradients existing at 1500-2500 m depth in the North Atlantic. Intermediate water is isotopically lighter in the South Atlantic than the North Atlantic, but these converge to approximately 0‰ in bottom waters (c.f. Bickert and Mackensen, 2003;Curry and Oppo, 2005). The very light values obtained in the Atlantic sector of the Southern Ocean may be influenced by a phytodetritus error (Mackensen et al., 1993), reflected in the large error bars. The data are not inconsistent with a decrease of δ 13 C with depth in the Pacific and Indian Oceans comparable to that observed in the South Atlantic, but there are insufficient data to establish a statistically significant vertical trend in these basins. Absolute δ 13 C derived from planktonic foraminifera is plotted for completeness, but these are subject to corrections for temperature and seawater chemistry. In Sect. 3.1.3, we described how Earth system model output fields could be used to correct for these, whereas proxy data are not generally available over 150 kyr. Therefore, these adjustments are not included in Fig. 10, and although there is reasonable spatial coherence in the Pacific and Indian basins, we do not use absolute planktonic δ 13 C in our interpretation. The deep ocean is isotopically heavier during the 5 ka time-slice (Fig. 11), indicating a decrease in the vertical δ 13 C gradient. Below 2500 m, mean Holocene LGM δ 13 C is 0.44 ± 0.13‰ in Indo-Pacific Ocean records, 0.57 ± 0.10‰ in North Atlantic records and 0.60 ± 0.18‰ in the South Atlantic records (errors are double the standard error, approximately 95% confidence intervals; see also Fig. 7). It is cautioned that this calculation is based on the assumption that errors are random, and that quoted errors apply only to locations where observations exist. Uncertainties in mean LGM δ 13 C within each ocean basin are much greater due to sampling biases. Nevertheless, the coherent signal observed throughout the deep ocean supports the hypothesis that the decrease in the vertical δ 13 C gradient is principally a global signal. At intermediate depths, LGM δ 13 C is typically indistinguishable from zero; however greater values of LGM δ 13 C are obtained in the East Atlantic than the West Atlantic.
Although most individual LGM δ 13 C observations cannot be distinguished from zero in the 47 ka time-slice (Fig. 12), the result that seawater was isotopically heavier during MIS 3 than during the LGM is statistically significant (mean LGM δ 13 C of 0.20 ± 0.06‰ in Atlantic records and 0.14 ± 0.11‰ in Indo-Pacific records). It is unclear whether this is a global signal, as four records suggest the presence of isotopically lighter waters in the intermediate-depth western Atlantic. There are no significant basin-scale differences between MIS 4 and LGM δ 13 C distributions (Supplementary Fig. 1). However, the greatest decrease in ocean δ 13 C occurred between MIS 5a and MIS 4. Below 2500 m in the 75 ka time-slice (Fig. 13), LGM δ 13 C is 0.20 ± 0.13‰ in Indo-Pacific records, 0.40 ± 0.08‰ in North Atlantic records and 0.36 ± 0.11‰ in South Atlantic records. The last interglacial (MIS 5e; 123 ka; Fig. 14) is indistinguish-able from MIS 5a in the Indo-Pacific (mean LGM δ 13 C of 0.17 ± 0.11‰) or the North Atlantic (0.44 ± 0.11‰), but there is evidence that South Atlantic waters were isotopically heavier (mean LGM δ 13 C of 0.53 ± 0.11‰) during the last interglacial than during MIS 5a, especially in bottom waters and Southern Ocean waters. Finally, an MIS 6 time-slice is provided in Supplementary Fig. 2, distinguishable from the LGM only by globally lower LGM δ 13 C values (mean LGM δ 13 C of −0.14 ± 0.10‰ in the Indo-Pacific, −0.15 ± 0.09‰ in the North Atlantic and −0.06 ± 0.09‰ in the South Atlantic).

Sensitivity to assumptions
In this section we examine the sensitivity of the data synthesis to the age model used and to assumptions made during the processing. All accompanying figures are provided in Supplementary Materials. Rather than propagating age error to an error in δ 13 C, we have quoted an estimated 2σ error of 6-8 ka on output ages. To investigate the effect of this on LGM δ 13 C, we randomly add noise, sampled from the distributions plotted in Fig. 4, to each record. Supplementary Fig. 3 shows the effect of this noise on mean and standard deviation of LGM δ 13 C timeseries in each ocean basin. Changes are typically less than 0.1‰, but can be greater for regions with few records such as the Indian Ocean. The addition of such noise does not significantly reduce the coherence of the global or regional signals. We also investigated the effect of regionally tuning age models to nearby 14 C-derived age models. This was done by calculating the mean offset between our age estimate and previously published 14 C-derived age estimates in each ocean basin, and adjusting the age of all records in that basin to compensate for this. These adjustments were only made for the period 0-50 ka, and were less than 4 ka at all time-slices in all regions. Supplementary Fig. 4 shows that the effect on LGM δ 13 C is negligible. We conclude that, because age errors are typically small compared with the timescale of interest, our findings are robust to the choice of age model. We did not apply a correction to LGM δ 13 C due to a time-variable phytodetritus effect, but instead added a random of error to records suspected of being influenced by this effect (flagged in Table 1). This leads to the hypothesis that the synthesis may be biased towards low δ 13 C at the LGM (high LGM δ 13 C at other times) by such records. In Supplementary Fig. 5, we compare the mean and standard deviation of the entire synthesis with records from upwelling regions, and find that the amplitude and phase of variability is similar. A major exception is the Holocene LGM δ 13 C in the South Atlantic, which is approximately twice as great in phytodetritus-flagged records than in the entire synthesis, caused by five records from the Atlantic sector of the Antarctic Circumpolar Current. There is no similar amplification of LGM δ 13 C variability prior to the LGM, suggesting that any isotopically light microenvironment that formed at these sites persisted for most of a glacial cycle. Attributing a greater uncertainty to these records, as we have done, is preferable to removing them from the synthesis.
Finally, we test the hypothesis that the amplitude of the glacial cycle is underestimated in records with low sediment accumulation rates due to filtering of variability by bioturbation in sediment. Comparing the mean and standard deviation of the entire synthesis with records from sites with accumulation rates less than 3 cm/ka ( Supplementary  Fig. 6), we find no evidence for lower amplitude variability in records from sites with low accumulation rates; in intermediate depths, the amplitude of the MIS 6 to MIS 5e deglaciation is typically greater in these records. It is likely that filtering due to bioturbation is most important on timescales that have already been removed by the application of a 10 ka filter. Therefore, the synthesis already excludes any millennial timescale peak in δ 13 C during interglacials, and this is not significantly exacerbated by the inclusion of records with low accumulation rates.

Discussion
We have presented benthic and planktonic δ 13 C data over the last 150 ka. The reasonable global coverage within the synthesis depends on the inclusion of records from a variety of species, from high-productivity regions, and from regions with low sedimentation rates. This approach makes uncertainties an essential part of the data set; these were obtained using the entire dataset to estimate the additional error present in data from low resolution records and in less reliable species. Errors are smaller for estimates of LGM δ 13 C, the anomaly relative to the local LGM value, than for the absolute δ 13 C because a large component of error is uniform within each record.
LGM δ 13 C time-slices have the advantage over δ 13 C time-slices as a modelling target that the latter target may favour states with errors that compensate for model error in the Holocene δ 13 C distribution. We caution that, due to imprecision in the age-models, the data synthesis should not be used to examine inter-basin differences on timescales shorter than 10 ka. A final caveat is that much remains to be learnt about how changes in properties other than seawater δ 13 C have influences the δ 13 C record. Several such properties (notably temperature, pCO 2 and carbonate ion concentration) are likely to change over a glacial cycle and introduce systematic errors to δ 13 C change estimates (Sect. 3.1.3). This is of particular importance for data from planktonic foraminifera. We therefore caution that recorded variability in planktonic δ 13 C may result from changes in ocean temperature and chemistry rather than in seawater δ 13 C.
When comparing the LGM time-slice with previous studies, it should be noted that a uniform adjustment, representing estimated changes in biosphere carbon storage (Duplessy et al., 1988b), is frequently applied to δ 13 C reconstructions, but is not applied in this study. We observe that isotopically light (δ 13 C < 0.5‰) bottom waters penetrate north as far as the Iceland Basin, attributed to an expanded Antarctic overturning cell in the Atlantic Ocean relative to today (Bickert and Mackensen, 2003;Curry and Oppo, 2005). These studies focused on compiling benthic LGM time-slices for the Atlantic, and included cores, mostly from the West Atlantic, that are absent from our synthesis due to the age modelling constraint (Sect. 2.1). The addition of further West Atlantic cores reveals a zonal gradient in the LGM South Atlantic at intermediate depths, with isotopically lighter water of southern origin in the west (Curry and Oppo, 2005), and isotopically heavier water of possible Mediterranean origin near the eastern boundary (Bickert and Mackensen, 2003). Our time-slice also excludes data from the Matsumoto and Lynch-Stieglitz (1999) compilation of southeast Pacific LGM data, where individual foraminer shells were interpreted as representing the LGM, based on δ 18 O measurements. They obtained glacial δ 13 C estimates of between −0.35 and −0.6‰ in the 2800-3800 m depth range. The result that the MIS 5e-MIS 6 δ 13 C difference is analogous to the Holocene-LGM difference supports the findings of Duplessy and Shackleton (1985), who described an increase in δ 13 C during deglaciation, although δ 13 C was generally lower during MIS 5e compared with the Holocene (Duplessy et al., 1984). Govin et al. (2009) examined changes in δ 13 C in 10 records from Atlantic and Antarctic source waters during glaciation. They found a two stage process, the decrease in Southern Ocean δ 13 C occurring early during glacial inception, and predating changes in the North Atlantic. We find this result to be robust, and that moreover there is negligible net decrease in δ 13 C outside the South Atlantic between the last interglacial and MIS 5a.
The data synthesis reveals a high degree of spatial coherence in δ 13 C variability in the deep Atlantic and Pacific Oceans, with high values during temperature maxima and low values during temperature minima. The amplitude of variability in the Atlantic Ocean is approximately double that in the Pacific, suggesting that the δ 13 C reservoir effect and ocean processes influencing deep Atlantic δ 13 C act in phase and are of similar magnitude. The relative lack of variability in intermediate-depth North Atlantic δ 13 C is consistent with a shoaling of 13 C-rich North Atlantic source waters during cold periods. The available data from outside the Atlantic Ocean also show less variability at intermediate depths. Atlantic δ 13 C appears to be influenced by increases in δ 13 C of the southern end-member during interglacials, but not during other temperature maxima. There is little evidence for basin-scale differences in δ 13 C between MIS 1, 5a, 5c and 5e, or between MIS 2, 4 and 6 (except for the globally lower δ 13 C values obtained during MIS 6), although individual records may show such differences. δ 13 C at a given location may be interpreted in terms of highly local processes rather than large scale processes discussed in this paper. We therefore recommend that care be taken when using small numbers of records to infer changes in global phenomena. The similarity between MIS 5a and 5e suggests that the relationship between δ 13 C and atmospheric pCO 2 over glacial cycles is complex. The major change in the δ 13 C record occurs between MIS 5a and MIS 4, whereas much of the decrease in pCO 2 (and almost all of the decrease in Antarctic temperature) occurs between MIS 5e and MIS 5a (Fig. 7). Therefore, either the drawdown of pCO 2 during early glaciation occurred by a mechanism that does not decrease deep ocean δ 13 C, or there was a compensating mechanism acting to increase deep ocean δ 13 C without increasing pCO 2 .
Constructing δ 13 C inventories would be a valuable tool in estimating changes in carbon storage in the terrestrial biosphere and shelf sediments. Duplessy et al. (1988b) used Pacific data only to estimate the increase in the mean δ 13 C of DIC in seawater between the LGM and the Holocene as +0.32‰. This yields a smaller change in carbon storage than is suggested by paleoecological reconstructions (Crowley, 1995;Peng et al., 1998), although the discrepancy would be partly explained if there was a lower proportion of isotopically light C 3 plants at the LGM than during the Holocene (François et al., 1999;Crucifix et al., 2005). Despite the coherence of the data presented here, we consider the coverage too incomplete to directly construct a time-series of δ 13 C inventories. The principal gaps in benthic data suited to our age modelling method are in the Indian Ocean, intermediate depths in the Pacific Ocean, the Southern Ocean, and pycnocline depths throughout the global ocean. Nevertheless, there are sufficient data to tightly constrain the evolution of an Earth system model, and we propose data assimilation into such a model as a viable means to reconstruct biosphere carbon storage over the last glacial cycle.