2.1 Carbonate burial in oceanic sediment

2.1.2 Production, export, and burial of carbonates in open ocean settings

In open ocean sediments, coccolithophorid and foraminiferal fragments are the main source of CaCO₃. The burial of CaCO₃ in deep-sea sediment is controlled by the production and export of CaCO₃ minerals from the surface ocean on one hand, and by the preservation of these phases at the seafloor on the other (Dunne et al., 2012). The CaCO₃ saturation state increases with the in situ concentration of Ca²⁺ and [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$], and decreases with increasing pressure and therefore water depth. The accumulation of respired carbon in deep water along the pathway of the global ocean circulation tends to acidify more poorly ventilated parts of the deep ocean, thereby reducing CaCO₃ preservation with water mass aging as shown in Fig. 2. Thus, CaCO₃ burial in the modern ocean is high in the Atlantic, where the carbonate compensation depth (CCD) is relatively deep, in the eastern equatorial Pacific, where production is high, and in the southern Indian Ocean, where high productivity coincides with elevated topographic features.

Figure 2(a) Deep-sea CaCO₃ burial based on sedimentary data; (b) CaCO₃ burial based on a metamodel optimized to fit the sedimentary data (Dunne et al., 2012). Filled circles show Holocene CaCO₃ burial from sediment cores used in the reconstruction of past burial fluxes, using the same color scale. Note that many sediment cores were collected at sites with accumulation rates that exceed the local average sedimentation rate, a common practice to obtain high-temporal-resolution records for paleoceanographic studies (see Sects. A1 and A4 in the Appendix).

Download

By contrast, the deep north Pacific, bathed by aged, DIC-rich waters, shows very low CaCO₃ burial rates. The fraction of pelagic CaCO₃ exported from the surface ocean that is eventually buried in sediments is estimated to vary between 0.6 and 0.06, depending on surface export, bottom water chemistry, sedimentary respiration rates, and lithogenic flux (Dunne et al., 2012).

Based on the modern surface sediment composition and estimates of bulk sediment burial rate, and also considering the export flux of CaCO₃ from the surface ocean, the modern CaCO₃ burial in deep-sea environments has been estimated to be between 100 and 150 PgC kyr⁻¹ (Catubig et al., 1998; Sarmiento and Gruber, 2006; Dunne et al., 2012). Thus, modern estimates suggest that carbonate burial is roughly equally distributed between shallow- and deep-sea environments but with significant uncertainty in their ratio (at least a factor of 2).

2.2 Organic carbon burial in marine sediment

2.3 Geological CO₂ emissions

Modern estimates of volcanic emission of CO₂ range between 6 and 14 PgC kyr⁻¹ for mid-oceanic ridges alone (Chavrit et al., 2014) and 44 PgC kyr⁻¹ for mid-oceanic ridges plus volcanic arcs (Coltice et al., 2004). Estimated emissions on geological timescales are notably higher, ranging from 40 to 150 PgC kyr⁻¹ (Sundquist and Visser, 2003; Burton et al., 2013), consistent with relatively quiet volcanic activity during the late Holocene (Gerlach, 2011). It has been suggested that increased emission during deglaciations, due to mantle decompression caused by ablation of glaciers and ice caps, could have reached 200 PgC kyr⁻¹ (Huybers and Langmuir, 2009; Roth and Joos, 2012). Continental rifting zones could also represent a significant source of exogenic carbon to the climatic system, amounting to 14 to 26 PgC kyr⁻¹ (Lee et al., 2016). Notably, a recent study significantly increased the potential contribution of solid Earth outgassing as compared to previous estimates (see details in Burton et al., 2013 and references therein), a revision that highlights the large degree of uncertainty in this important flux.

Some reconstructions for past global volcanic activity have suggested a tight coupling between glacial ice volume and volcanic eruption frequency (Jellinek et al., 2004). Increased eruption frequency among subaerial volcanoes appears to have occurred following the highest rate of sea level rise (with a 4 kyr lag), presumably in response to mass unloading and isostatic adjustment (Huybers and Langmuir, 2009; Roth and Joos, 2012; Kutterolf et al., 2013). Subsea volcanism at mid-ocean ridges also appears to vary in response to isostatic loading as ice sheets grow and sea level drops. Lund et al. (2016) suggested that the maximum rates of seafloor CO₂ emission could occur approximately 10 kyr after a sharp sea level drop, which occurred most recently at the onset of the Last Glacial Maximum (LGM) (Lambeck and Chappell, 2001; Waelbroeck et al., 2002; Grant et al., 2012), implying that peak seafloor emissions would be roughly coincident with the start of deglaciation. The isotopic composition of volcanic emissions has been estimated at around −5 ‰ (Kump and Arthur, 1999; Cartigny et al., 2001; Deines, 2002; Roth and Joos, 2012).

2.4 Erosion of continental surfaces

The second main source of geological carbon to the active pool consists in the transfer of carbon from exposed rocks through physical and chemical weathering. Carbonate rock weathering has been estimated to currently add between 70 PgC kyr⁻¹ (Hartmann et al., 2014) and 155 PgC kyr⁻¹ (Gaillardet et al., 1999) to the active pool, while the petrogenic particulate C_org flux to the ocean is thought to amount to between 40 PgC kyr⁻¹ (Copard et al., 2007; Galy et al., 2015) and 70 PgC kyr⁻¹ (Petsch, 2014). Alkalinity fluxes from the continent to the ocean also impact the carbon system. Most of the alkalinity released from the continental surface corresponds to the weathering of carbonate (63 to 74 %; Gaillardet et al., 1999), with the remainder supplied by silicate weathering. Total alkalinity fluxes to the ocean amount to between 30 Pmol equ kyr⁻¹ (Amiotte Suchet et al., 2003) and 40 Pmol equ kyr⁻¹ (Gaillardet et al., 1999).

The degree to which weathering varied over glacial–interglacial cycles is highly uncertain. An early modeling study estimated that global chemical erosion was 20 % higher during the LGM (Gibbs and Kump, 1994), mostly because of enhanced erosion of newly exposed continental shelves. However, Vance et al. (2009) have argued that silicate chemical weathering was much lower during glacials. It is commonly assumed that physical erosion of continental surfaces was more rapid during glacials because of enhanced action of ice sheets at high latitude (Herman et al., 2013) and the exposure of shelves, although this too has been debated (von Blanckenburg et al., 2015). It is also possible that chemical weathering could have been reduced, or changed in nature (Torres et al., 2017), due to lower temperature so that global erosion remained relatively stable over glacial–interglacial timescales (Foster and Vance, 2006).

The isotopic composition of the global erosion flux is mostly determined by the ratio of C_org to carbonate weathering and has been estimated to range between −8.5 ‰ and −5 ‰ (Kump and Arthur, 1999; Roth and Joos, 2012).

3.1 Analytical strategy overview

Sediment cores sample only a tiny fraction of the global seafloor. In order to use these unevenly distributed observations to infer integrated changes in global carbon burial, we need to take into account the following aspects: (1) the limited spatial extent that is integrated by each core; (2) the fact that accumulation rates at an individual site may not be representative of the accumulation within a broader region; and (3) the fact that carbonate accumulation rates vary widely between different parts of the world.

To account for these potential shortcomings, our quantification of global carbonate burial includes three elements:

• The first element is an estimate of the spatial variations in global carbonate mass accumulation rate (MAR) in the modern ocean derived from sediment core tops and process modeling (Dunne et al., 2012).

• The second element is subdivision of the ocean into coherent biogeochemical provinces, in order to account for the spatial and bathymetric heterogeneities in the sedimentary carbonate burial record. Because they take into account spatial structure in carbonate production and preservation, the province approach is likely to be an improvement over simple interpolation, especially in regions where measurements are scarce (see details in Appendix A), though future work may be useful to quantify exactly how the approaches differ.

• The third element is reconstruction of the relative variation of the carbonate MAR in each province over time, based on a compilation of available marine sediment proxy records. These relative changes are then multiplied by the modern CaCO₃ MAR estimate in each province and integrated to assess the global history of deep-sea carbonate burial (see also Cartapanis et al., 2016 for further details). We note that the last step implicitly assumes that carbonate currently present in modern sediment is in equilibrium with its surrounding environment and is not being gradually dissolved (Berelson et al., 2007) or precipitated (Sun and Turchyn, 2014), which is not universally correct. In the future, this limitation could perhaps be overcome by assimilating the observed MAR into a mechanistically based sediment module coupled with carbonate production and ocean chemistry, though this is beyond the scope of the present work.

Figure 3Coverage and number of records used for each province map. Coverage is calculated by summing the modern burial rates in each province (Fig. A1) for which at least one sedimentary record is available, divided by the modern burial flux. The modern map used corresponds to Fig. 2a.

Download

3.2 Modern carbonate MAR map

We use two different but closely related maps of modern carbonate MAR (Fig. 2), both of which were previously published by Dunne et al. (2012). The first is based only on observations and was determined by multiplying the carbonate concentration in surface sediments by the bulk sediment accumulation rate. The carbonate concentration map was generated by combining several literature sources (Seiter et al., 2004; Dunne et al., 2007), whereas the bulk MAR map was determined based on the geometric mean of two pre-existing maps (Jahnke, 1996) and other unpublished data; see details in Dunne et al. (2007, 2012). It is important to note that the observationally based estimates of both the carbonate concentration and MAR have insufficient spatial resolution to adequately resolve continental margins (Fig. 2a). The second map was generated using a metamodel derived from a process-based approach (including satellite-derived carbonate export production and water-chemistry-dependent dissolution in the water column and sediments) that was minimally tuned using the sedimentary data (Dunne et al., 2012). The metamodel map shows similar patterns in the deep sea but better resolves the bathymetric imprint on carbonate preservation (Fig. 2b). The two approaches provide slightly different spatial distributions of burial, allowing a test of the sensitivity of our reconstructions to this uncertainty (see Table 1). Although the absolute values differ, the relative changes in our reconstructed global carbonate burial over time are insensitive to the choice of maps (see Appendix A).

Table 1Estimates of the modern global deep-sea CaCO₃ burial flux. The fractional contribution of each basin is shown, in percent, and the total estimated burial flux shown in g CaCO₃ yr⁻¹. The two columns correspond to the two estimates shown in Fig. 2, and the basin areas correspond to province map 1 shown in Fig. A1.

Download Print Version | Download XLSX

The two modern maps attempt to accurately represent the spatial distribution of the modern CaCO₃ burial flux. Sediment cores, however, are often situated to target rapidly accumulating depocenters in order to provide greater temporal resolution for paleoceanographic reconstructions. As a result, we expect a bias towards high accumulation rates, which is confirmed in Fig. 2. We address this sampling bias when evaluating uncertainties (see Sects. A1 to A4).

Figure 4Example of calculated changes in deep-sea sedimentation relative to the Holocene using one province map. The first three rows show relative changes in bulk sediment MAR, carbonate content, and carbonate MAR, while the last row shows absolute changes in carbonate MAR. Changes are shown for the LGM and Marine Isotope Stages 4, 5e, and 6. Note that the absolute carbonate flux refers to the integrated burial over the province (not area normalized). The modern map used here corresponds to Fig. 2a, and the province map is S.L.2. Note that the data shown here are not corrected for the long-term trend of accumulation rate, as the correction is calculated and applied at the global burial level.

Download

3.3 Relative changes in burial rate from sediment core database

We created a database including surface- and downcore sediment datasets, retrieved from the National Oceanic and Atmospheric Administration (NOAA) (ftp://ftp.ncdc.noaa.gov/pub/data/paleo/paleocean/sediment_files/, last access: Mai 2015) and PANGAEA (http://www.pangaea.de, last access: Mai 2015) online repositories. All the datasets that contained CaCO₃ concentrations, CaCO₃ accumulation rates, age models and sediment density values were taken into consideration (for a complete detailed description of the procedure and quality control, the reader is referred to Cartapanis et al., 2016). Whenever available, ²³⁰Th-normalized sediment accumulation rates were used. Where these measurements were not available, sedimentation rates were determined at the depth of each carbonate measurement from the reported age–depth relationship. When possible, we used reported dry bulk density values. Otherwise, we assumed a constant value of 0.9 g cm⁻³ to convert linear accumulation rates to MAR. Each sedimentary record was expressed as a fraction of the respective mean Holocene value.

The mean LGM : Holocene (18–25 and 0–10 ka, respectively) ratio was determined for each province and multiplied by the corresponding modern burial rate to estimate the absolute value of carbonate burial during the LGM. Note that we used the geometric mean of available records in each province, given that the LGM : Holocene ratios of individual cores are log-normally distributed. Moreover, if in a single province we have two MAR records, one showing 2 times higher MAR during the glacial, while the other showed only half of modern during the glacial, the best guess would be that the true change in sedimentation was negligible. This is better represented by the geometric mean of (2, 0.5) = 1 (no change), whereas the arithmetic mean (2, 0.5) = 1.25 would bias the result toward high values. The carbonate burial of all provinces is summed to provide the LGM global burial rate.

The same procedure was applied to calculate global carbonate sediment content and bulk MAR at every 1 kyr time step over the past 150 kyr. Any time interval within a sediment record for which no measurement was present over a period longer than 10 kyr was excluded from the calculation of the regional stacks, leading to a total of 637 sediment time series, from which 370 records were used to calculate the LGM : Holocene ratio (Figs. 3–4), while 170 to 130 records were used for periods older than 40 kyr. When no record was available for a specific province and/or for a specific time period, the province flux was assumed to be constant and equal to modern values. The provinces for which no downcore record was available represent 0 % to 18 % of the ocean surface, corresponding to 0 % to 12 % of modern burial (this metric will be referred as “coverage” in the following discussion) for the LGM (18–25 kyr; Fig. 3). As a consequence, the reduction of the coverage during older periods tends to diminish the difference from the modern situation, a conservative assumption that may bias our reconstruction toward smaller changes than actually occurred. The available records are distributed over the entire depth range (100 records shallower than 1000 m, 103 records between 1000 and 2000 m, 137 records between 2000 and 3000 m, 290 records between 3000 and 4000 m, and 106 records below 4000 m deep).

Figure 5CaCO₃ MAR time series, considering different aspects of variability. (a) Bulk sediment MAR time series distribution based on the 20 province-map scenarios (mean, black line; ±1σ, red area; ±2σ, blue area), and the distribution of the long-term trend correction for each province map (mean, black line; ±1σ, dark grey; ±2σ, light grey) using a least-squares spline model. (b) Bulk sediment MAR as in panel (a) but with the modeled long-term trend removed. (c) Global CaCO₃ MAR reconstruction without detrending. (d) Global CaCO₃ MAR reconstruction calculated only from variations in CaCO₃ concentration, assuming a constant bulk sedimentation rate. (e) Global CaCO₃ MAR reconstruction obtained from the bulk MAR and variations in CaCO₃ concentrations. Panel (f) is the same as (e) but with MAR corrected for the long-term trend. In all panels, the thin black lines correspond to the 20 individual province-map scenarios. The modern map used for all calculations corresponds to Fig. 2a.

Download

Figure 6Principal components of temporal variation in sediment mass accumulation rates and CaCO₃ concentrations. The first principal component (PC1) is displayed in panel (a) for both sedimentation rates and carbonate content, and explains 28.6 % and 33.1 % of the total variance, respectively. Yellow vertical bars correspond to Holocene and MIS5e; blue vertical bars are for LGM, MIS4, and MIS6. Panels (b) and (c)show the correlation coefficients for individual records of sedimentation rates and carbonate content, against their respective PC1. In this analysis, we used only records that cover 80 % of the studied interval.

Download

3.4 Assessing changes in deep ocean MARs and their role on the global estimates

Since carbonate MAR depends on both the sedimentation rate and sedimentary carbonate concentration, we first calculate the global bulk MAR by considering only relative changes in the sedimentation rates and a modern bulk MAR map (Jahnke, 1996). This reconstruction (Figs. 4, 5A, and Table 3) shows generally higher bulk sediment burial during peak glacials and deglaciations, broadly consistent with previous studies (Covault and Graham, 2010). However, the reconstructions display a pronounced increasing trend towards the present for all the province-map scenarios (Fig. 5A). This is partly due to the absence of density data for many cores, which tend to underestimate burial for older, more compacted sediment. Moreover, the apparent sedimentation rate decreases as a power–law function of the intervals between two age controls (Schumer and Jerolmack, 2009), suggesting that sedimentation rates are likely underestimated as age increases, when the age model accuracy is generally lower compared to more recent periods.

Because of the lack of density information for individual cores, we chose to correct for compaction at the global scale by assuming that the global bulk burial rate was equal during the late Holocene and Marine Isotope Stage 5e (MIS5e), and subtracting the intervening trend, an assumption that could be revisited in future work as more information becomes available. Given this assumption, we corrected the long-term trend of the global aggregate bulk sedimentation rate for each province map using a least-squares spline (Fig. 5a and b). The detrended CaCO₃ MAR reconstruction (Fig. 5f) is much more similar to the globally flux-weighted CaCO₃ concentration record (Fig. 5d) and suggests that global bulk MAR in deep-sea sediment peaked during glacial maxima and deglaciations, as well as during MIS4: 114 %, 104 %, and 105 % relative to the Holocene MAR for the LGM, MIS4, and MIS5e (119–124 ka), respectively. The LGM-to-Holocene change is consistent with, but smaller than, a recent estimate (125 % ± 15 %) using a smaller but partially overlapping set of cores that included C_org measurements (Cartapanis et al., 2016).

We performed a principal component analysis (PCA) on the MAR time series in order to assess the overall spatial and temporal variability (Fig. 6). For this calculation, we used only the records that cover at least 80 % of the studied interval. The first principal component (PC) of the bulk MAR shows variability that is similar to the bulk MAR reconstruction, and the sedimentation rate time series are broadly correlated with the sedimentation rate PC1, without showing any consistent geographic pattern.

We then performed a similar analysis, using only the relative changes in the sedimentary carbonate concentration (assuming constant sedimentation rates), and the modern carbonate MAR map as a reference (Dunne et al., 2012). This analysis shows that the CaCO₃ content in deep-sea sediment (weighted by modern carbonate MAR) was consistently lower during glacials compared to interglacials, with the lowest values reported during MIS4 and MIS2 (Figs. 4, 5d, and 6). The first PC of the carbonate content records is similar to the global reconstruction, and the correlation map between CaCO₃ concentration time series and CaCO₃ concentration PC1 shows a clear antiphasing between the eastern Pacific and the rest of the world ocean (Fig. 6). While CaCO₃ concentrations were lower during the glacial in the Atlantic, they were substantially higher in the eastern Pacific.

In order to assess the impacts of variations of both bulk MAR and CaCO₃ concentration in our reconstruction, we calculated global carbonate MAR variations using two different methods. First, we calculated global CaCO₃ MARs using the relative changes in individual sedimentary carbonate MAR time series (Figs. 4, 5c). Second, we made the simpler calculation of taking the product of the reconstructed CaCO₃ concentration with the average relative changes in bulk MAR reconstruction for each scenario (Fig. 5e).

As both methods provided very similar results, we calculated CaCO₃ MAR using the simpler of the two, taking the product of the relative changes in corrected bulk sediment MAR (Fig. 5b) and the reconstruction of CaCO₃ concentration (Fig. 5d). This final reconstruction (Fig. 5f; see details in Fig. A4) thus accounts for variations in the bulk MAR, as well as variations in the concentration, and is corrected for the long-term trend directly inferred from the bulk MAR.

As shown in Fig. 5f, the reconstruction indicates that global mean deep-sea CaCO₃ burial was not significantly different from the present day during MIS5, dropped to the lowest rates during MIS4, and then remained below interglacial rates until the deglaciation. In part, the similarity of CaCO₃ burial between MIS5 and the Holocene reflects our assumption that bulk sediment MARs during MIS5 and the Holocene were similar, which underlies our detrending. However, MIS5 and the Holocene are also indistinguishable in the burial flux reconstruction estimated from CaCO₃ concentration alone (Fig. 5d), in which MIS4–MIS2 also stand out as having significantly reduced burial.

4.1 Mass-balance calculation

In order to estimate the impact of the reconstructed carbonate flux on the carbon and alkalinity budgets, we made a simple running-sum calculation of carbon and alkalinity sources and sinks over the full last glacial cycle. We did so within the numerical framework of a simple three-box model (three oceanic boxes and the atmosphere) using the input and output fluxes for both carbon and alkalinity as illustrated in Fig. 1 and model settings described in Toggweiler (1999) following the equations detailed in Sarmiento and Toggweiler (1984). Briefly, the model tracks the mass budget for inorganic carbon in three well-mixed ocean boxes (surface high and low latitudes, and deep ocean) and in the atmosphere. The C cycle is linked to the cycles of alkalinity and nutrients (phosphate) in the ocean boxes. For simplicity, nutrients are prescribed in the surface boxes, and water mass transport between boxes is kept constant. Exchange of CO₂ between the surface boxes and the atmosphere exchange is based on aqueous carbon chemistry and a constant piston velocity (Sarmiento and Gruber, 2006). The deep-sea [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] concentration was calculated based on aqueous carbon chemistry (Sarmiento and Toggweiler, 1984; Dickson and Goyet, 1994) and can be very well approximated by the difference between dissolved inorganic carbon and alkalinity. We emphasize that the use of the three-box model here only provides a very small modification of the results obtained with the simple running sums, providing average deep ocean values that vary slightly due to changes in air–sea partitioning; very similar results would be obtained by considering simply the mass balance of global alkalinity and carbon inventories, without calculating their distributions between ocean and atmosphere.

The two input fluxes to the active pool are the geological outgassing of CO₂ and the erosion of continental surfaces (due to the dissolution of carbonate and oxidation of sedimentary C_org), while the four output terms represent carbon burial in sediment divided into four components: both C_org and CaCO₃ for each of the shelf and deep-sea environments. The input flux of alkalinity is assumed to be due to the erosion of continental surfaces, including both carbonate and silicate weathering, while the flux out of the system is given by burial of CaCO₃ in marine sediment.

Our mass-balance simulations prescribe a range of possible scenarios for variable input and output fluxes, in order to estimate the resulting transient changes of carbon and alkalinity inventories and δ¹³C. The deep-sea burial fluxes are derived from our sedimentary reconstructions, while the shelf burial fluxes and the input flux are varied in order to illustrate hypothetical (yet realistic) scenarios. In total, we include eight different scenarios to illustrate the impact of each variable in isolation, as well as in combination with others (Table 3). In all cases, we vary the fluxes over feasible ranges using a Monte Carlo approach (1000 transient simulations per experiment).

It is important to note that we do not consider temporally variable erosion fluxes but instead hold them constant over time in any given simulation, due to a paucity of consistent observational constraints. The work of von Blanckenburg et al. (2015) argued that erosion was, indeed, constant over glacial cycles, while others have argued that erosion was lower (Vance et al., 2009) or higher (Herman et al., 2013) during glacials. Erosion contributes both carbon and alkalinity, and although it remains difficult to evaluate the global impacts of changes in erosion, they may have been substantial.

A key feature of the simulations is that they assume long-term balance over a full glacial cycle. The fact that CO₂ and δ¹³C display very stationary mean values over several glacial cycles (Zeebe and Caldeira, 2008) suggests that the system is indeed close to steady state over each full glacial cycle, so that the assumption is unlikely to introduce a large bias. Therefore, assuming that the system was at steady state over the last glacial cycle, we calculated the input fluxes of carbon and alkalinity that would have been necessary to balance the prescribed burial (as well as prescribed input changes, if any) over the full glacial cycle for a given scenario and provided those as a constant input flux throughout the cycle. Thus, although the carbon and alkalinity fluxes are out of balance at each point within the cycle, they always integrate to zero net change between 0 and 123 kyr. Similarly, the isotopic composition of the system is prescribed by setting the δ¹³C of the erosion flux at a value that balances the other fluxes, consistent with the long-term steady-state assumption. Although forcing the carbon isotopes to be equal at start and end points does ignore a small long-term drift in whole-ocean δ¹³C (Oliver et al., 2010), it provides a clearer illustration of the isotopic implications of inventory changes over a single glacial cycle. We note that this strategy results in different total input fluxes (and therefore different residence times of the active pool) for each experiment.

To account for uncertainties in the modern fluxes (Fig. 1), the relative changes in the prescribed time series are multiplied by an estimate of the global modern flux, which is generated using a distribution that represents the uncertainty in modern global burial estimates (see Table 3, and Appendix B and Fig. A6). The forcing thus accounts for our reconstruction and their uncertainties, as well as the uncertainties in the modern burial estimates. The same procedure was applied for C_org burial variations in the deep ocean as shown in Cartapanis et al. (2016) but using a different distribution (Table 3; see details in Fig. A6).

In the following sections, we consider eight different scenarios in which we calculate the total mass and δ¹³C of the active carbon pool, as well as the oceanic alkalinity inventory, according to different combinations of temporally varying input and output fluxes for carbon and alkalinity. In experiments 1 and 2, we isolate the impact of reconstructed deep ocean inorganic (this study) and organic carbon burial (Cartapanis et al., 2016), respectively, and combine them in experiment 3. Experiment 4 only considers scenarios for shelf burial variations, while experiment 5 combines experiments 3 and 4, thus including variations of burial in both deep and shallow environments. Experiment 6 depicts a hypothetical deglacial carbon pulse based on Huybers and Langmuir (2009) and is combined with experiment 3 (burial in the deep sediment) in experiment 7, and with experiment 5 (burial in deep and shallow sediments) in experiment 8. Thus, these scenarios incrementally explore the impact of the different forcing and their interactions on the carbon and alkalinity budgets over the last glacial cycle.

Table 2Summary of estimated changes in carbonate MAR for each province map. The number of records used, mean coverage (%), and mean MAR (in % of modern value) are shown for four intervals (LGM, MIS4, MIS5e, and MIS6) for each province map. The results correspond to the calculations of temporal variability shown in Fig. 5f. The modern burial map used for all calculations corresponds to Fig. 2a.

Download Print Version | Download XLSX

Table 3Modern burial distribution (PgC kyr⁻¹) used for the Monte Carlo simulations (mean and standard deviation). When the distribution type is not specified, a normal distribution was used. Boldface indicates variable fluxes over time. Results of experiments 1 and 2, 3 to 5, and 6 to 8, are shown in Figs. 7, 8, and 10, respectively. The gamma distributions (γ distribution) used for C_org burial are shown in Fig. A6.

Download Print Version | Download XLSX

4.2 Deep-sea burial impact on carbon and alkalinity inventories (experiments 1 to 3)

Assuming a constant alkalinity flux to the ocean, and all other factors being equal, our reconstruction of deep-sea CaCO₃ burial in the deep ocean coupled with uncertainties in modern fluxes (Table 2) implies changes in the alkalinity inventory that could reach up to 100 µequ kg⁻¹, with a minimum centered around ∼70 kyr near the MIS4/MIS5 transition, gradually rising to a maximum during the early deglaciation at ∼15 kyr (Fig. 7a, experiment 1). This increasing trend is a robust consequence of the low reconstructed carbonate burial from MIS4 to MIS2, with a most probable amplitude of 70±30 µequ kg⁻¹. Note that the [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] changes shown in all of our experiments ignore possible changes in ocean circulation or biology, and they do not consider contributions to alkalinity other than ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$, and are therefore not expected to be accurate – they are simply illustrations of the purely inventory-driven expectations of deep-sea [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] change. The actual realized changes of [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] that occurred over the last glacial cycle have been shown to have been quite small (Yu et al., 2010), which would reflect the net result of any inventory changes plus other processes such as ocean circulation and biology-driven changes.

Changes in C_org burial in the deep sea followed a markedly different progression, with larger (rather than smaller) burial during peak glacial and deglaciation (Cartapanis et al., 2016), causing the C inventory to rise over MIS5 by 400 (150 to 950) PgC, equivalent to 20 (10 to 50) µmol kg⁻¹ in the global ocean (Fig. 7b, experiment 2), and to decrease by the same amount from MIS3 to the early Holocene. Summing the two deep-sea burial fluxes together leads to a slightly muted change in the carbon inventory (Fig. 8a, experiment 3). The changes in deep-sea C_org burial would have decreased the isotopic composition of the active carbon pool during the early glacial, by 0.2 (0.1 to 0.5) ‰ during MIS4, but with a negligible change between LGM and the late Holocene (Fig. 8a). It is important to note that, for these three first experiments, the mean C input flux required to compensate for the outputs is close to 500 PgC kyr⁻¹, which exceeds the standard estimates for present-day carbon input (160–400 PgC kyr⁻¹; Fig. 1). This discrepancy arises from the relatively large estimates of present-day carbon burial, combined with the absence of much lower deep-sea burial fluxes during the glacial – either the present-day input fluxes are underestimated, the present-day burial is overestimated, or there were additional changes in the fluxes over the glacial cycle. Thus, we turn next to possible changes in burial on the shelves.

Figure 7Simulated impacts of reconstructed deep burial fluxes on carbon and alkalinity inventories: (a) variable deep-sea CaCO₃ burial and (b) variable deep-sea C_org burial. For each panel, the uppermost curves show the carbon budget terms (constant background carbon input, organic, inorganic carbon burial, and net flux from the geological to the active carbon reservoir) used to force the simulations. For each time series, the mean of the scenarios (thick line), ±1σ, and ±2σ are shown. Below are shown the simulated δ¹³C of the total active carbon pool, deep ocean alkalinity, deep ocean DIC, and implied deep ocean [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] (µmol kg⁻¹). The [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] only reflects changes in the alkalinity and DIC inventories, not changes in ocean circulation or biology. An estimate for the whole-ocean δ¹³C by Peterson et al. (2014) during the LGM is also shown in orange (with upper and lower estimates). Vertical bars highlight cold periods (MIS2/LGM and MIS4) in blue and full interglacial conditions (Holocene and MIS5e) in yellow. The distribution of the modern fluxes used to drive the scenarios is shown by the histograms at the bottom of each panel.

Download

Figure 8Simulated impacts of deep and shelf burial fluxes on carbon and alkalinity inventories: (a) variable deep-sea burial (carbonate plus C_org), (b) variable shelf burial, and (c) variable total marine burial (deep plus shelf). All quantities and bars are plotted as in Fig. 7.

Download

4.3 Shelf burial impact on carbon and alkalinity inventories (experiments 4 and 5)

The burial of CaCO₃ and C_org on continental shelves depends both on their production rate and their preservation. Preservation is particularly difficult to estimate, given the importance of sediment transport and redeposition in dynamic coastal environments. Here, we use the very simple illustrative assumption that, globally, the burial over shelves was linearly proportional to the changes in the surface area available for sediment production and accumulation. Although this is a gross simplification of reality, and may overestimate the actual changes, it allows a first-order estimate of the potential magnitude of shelf burial changes on the global carbon cycle (Wallmann et al., 2016).

We calculated the relative changes in the surface area bounded by the 0 and 100 m isobaths as well as the 0 and 200 m isobaths over time, using ETOPO01 (Amante and Eakins, 2009) and a global sea level reconstruction curve (Grant et al., 2012). Because of the distribution of the global hypsometric curve (Fig. 9a), the surface area of the 0–100 m depth range was reduced by up to 60 % during the last glacial, which is certain to have led to reduced burial (Fig. 9b), consistent with Vecsei and Berger (2004).

Our results show that variations of C_org and CaCO₃ burial fluxes driven by changes in sea level could have increased DIC and alkalinity inventories by up to 150 µmol kg⁻¹ and 150 µequ kg⁻¹ from MIS5 to the last deglaciation, because of reduced burial during times of lower sea level (Fig. 8b, experiment 4). The ALK and DIC inventories change in parallel because of our assumption that both burial fluxes would have varied in proportion to the shelf area, which leads to relatively small changes in the [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] (the difference between DIC and ALK). Despite these parallel changes, the fact that C_org burial is strongly focused in shallow environments, in contrast to the equal partitioning of carbonate burial between shallow and deep environments, means that the global proportion of C_org to carbonate burial would have varied with changing shelf area (Wallmann et al., 2016). Our calculations suggest that the relative reduction of C_org burial cause by limited shelf area during glacials would have lowered the δ¹³C of the active carbon pool by up to 1 ‰ from MIS5 to the early deglaciation. The large uncertainties in modern burial fluxes result in large uncertainties in the ratio of DIC : ALK changes, with important implications.

Figure 9(a) Global hypsometric curve based on ETOPO01 (Amante and Eakins, 2009); (b) LR04 benthic stack (Lisiecki and Raymo, 2005); sea level reconstruction (Grant et al., 2012); and calculated changes in the surface of the 0–100 and 0–200 m depth intervals relative to the modern situation. Vertical bars are as in Fig. 6.

Download

4.4 Combined impact of sedimentary burial and carbon input variations (experiments 6 to 8)

It has been suggested that the rate of carbon input to the active pool from volcanic sources increased during glacial maxima (Lund et al., 2016; Ronge et al., 2016) and/or deglaciations (Huybers and Langmuir, 2009; Stott et al., 2009). Note that, for the purpose of this discussion, we consider peatlands and permafrost part of the “active” pool but not surficial marine sediments. The magnitude of these postulated volcanic carbon input variations is highly uncertain, and reconstructions of relative changes in subaerial volcanic outgassing rates over time suggest that the volcanic carbon input was notably delayed, relative to the deglacial increase in PCO₂ recorded in ice cores, requiring other mechanisms to be important at this time (Huybers and Langmuir, 2009; Roth and Joos, 2012). Nonetheless, it would appear likely that changes in geological carbon inputs over the glacial cycle were significant. We do not attempt to provide an exhaustive test of the existing hypotheses for temporally variable carbon input. Instead, we simply apply a pulse of carbon input to our mass-balance calculation during the deglaciation, which follows the temporal changes given by the reconstruction of changes in terrestrial volcanic emission (Huybers and Langmuir, 2009) to produce an illustrative variable carbon input to the active pool.

As shown in Fig. 10a (experiment 6), the simulation with only variable carbon input produces similar results over the deglaciation to those estimated by Huybers and Langmuir (2009) (see also Roth and Joos, 2012). However, within our broader perspective, it is clear that, when viewed over the glacial cycle (and given our requirement that the carbon inventory be the same at the beginning and end of the cycle), a major effect of a temporally variable carbon input is to lower the implied background input rate. As a result, the carbon inventory gradually declines from interglacial to glacial maximum, and is then rapidly increased by the volcanic pulse (Fig. 10a).

When our reconstructed deep burial variations are also included, the magnitude of the MIS5–MIS4 [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] peak is reduced relative to the deep-burial-only case (Fig. 10b, experiment 7). When shelf burial changes are also included, the mean carbonate ion variations are altered little, but the δ¹³C changes significantly (Fig. 10c, experiment 8). In this case, the δ¹³C of the ocean during the LGM is similar to that observed, even with no change in the terrestrial biosphere, as discussed further below.

Notably, only the experiments including variable shelf burial (Figs. 8b, c, and 10c) simulate Holocene input fluxes that overlap with the observational estimates summarized in Fig. 1. This is because only the variable shelf experiments provide a sufficient reduction of carbon burial fluxes during glacial periods to allow the full cycle carbon budget to be balanced with the relatively small modern estimates of the carbon input fluxes. Although this should not be taken as proof that the changes are correct, it shows encouraging consistency with the reconstructed and conjectured evolution of the carbon inventory over time.

Figure 10Simulated impacts of variable burial and outgassing fluxes on carbon and alkalinity inventories: (a) deglacial carbon input pulses, (b) variable deep-sea burial plus deglacial carbon input pulses, and (c) all variable (deep-sea burial, shelf burial, and deglacial carbon input pulses). All quantities and bars are plotted as in Fig. 7, with the addition of time-varying carbon input pulses based on subaerial volcanic emission (Huybers and Langmuir, 2009) in red.

Download

5.1 Carbonate burial in the deep ocean

The calculations shown above are based on individual time series of CaCO₃ MAR from hundreds of sediment cores, using previously reported sedimentary CaCO₃ content, sedimentation rates, and dry bulk density when available (see Supplement). Following the method of Cartapanis et al. (2016), we used these individual records to calculate, at each 1 kyr time step over the past 150 kyr, the CaCO₃ MAR in oceanographically defined provinces as a fraction of their Holocene values. Notably, these sediment cores record the actual realized burial, not just a proxy thereof. Although the preserved record may underrepresent transient deposition and dissolution events, and the available sediment cores are undoubtedly subject to multiple biases (Appendix A), our broad spatial reconstruction should tend to minimize locally obfuscating factors, providing a more robust reconstruction of carbonate burial, globally. Absolute changes in CaCO₃ burial were obtained by scaling reconstructed relative changes in CaCO₃ MAR to the modern MAR (Figs. 2 and A1). Uncertainty in the global burial flux was estimated using a bootstrap technique (Appendix A). The resulting reconstruction is sure to contain biases and will hopefully be improved in future; however, it represents the most comprehensive such undertaking to date and it is not obvious to us how the reconstruction could be improved significantly without a major expansion of the existing dataset.

As shown in Fig. 11, the reconstruction indicates that global mean deep-sea CaCO₃ burial was not significantly different from the present day to MIS5, dropped to the lowest rates during MIS4, and then remained below interglacial rates until the deglaciation. The spatial pattern of change, which displays the well-known deep Pacific–Atlantic antiphasing (Farrell and Prell, 1989; Hain et al., 2010) (Fig. 6), is consistent with a previously reported shoaling of well-ventilated North Atlantic Deep Water (NADW) (Boyle and Keigwin, 1982; Curry et al., 1988; Labeyrie et al., 1992; Lynch-Stieglitz et al., 2007; Guihou et al., 2011) and an enhanced soft tissue pump at the start of MIS4 (Kohfeld and Chase, 2017). The associated increase of DIC-rich Antarctic Bottom Water (AABW) would have decreased [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] and favored the dissolution of CaCO₃ in deep Atlantic sediments. An increased rain of organic matter to the deep sea, or a change in the rain ratio, such as what may have resulted from a postulated shift toward proportionally more siliceous phytoplankton (Brzezinski et al., 2002) or a midlatitude temperature decrease (Dunne et al., 2012), could have accentuated the DIC increase (Boyle, 1988; Hain et al., 2010; Yu et al., 2014; Galbraith and Jaccard, 2015).

Consequently, the reduced rate of alkalinity removal in the Atlantic during the glacial would have increased the total alkalinity inventory of the ocean, which would be expected to have increased [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] and consequently carbonate preservation elsewhere, a compensation long understood to have occurred in the Pacific Ocean (Farrell and Prell, 1989) and consistent with our observed Atlantic–Pacific antiphasing. However, our calculations suggest that reduced burial in the Atlantic Ocean during the last glacial was not entirely compensated in the Pacific Ocean, leaving glacial carbonate burial rates lower than for the Holocene (Figs. 6 and 11). This implies that, if the alkalinity inputs from weathering had remained constant (as assumed in our mass-balance calculations), excess alkalinity in the global ocean would have built up during glacials as compared to interglacials. Additionally, despite notable pulses of CaCO₃ burial in the deep ocean in some locations during the deglaciation (Jaccard et al., 2005; Dubois et al., 2010; Rickaby et al., 2010; Jaccard et al., 2013), and in contrast with expectations (Broecker and Peng, 1987; Boyle, 1988; Sigman and Haug, 2003), our reconstructions suggest that preserved deep-sea CaCO₃ burial itself did not drastically peak during the last deglaciation, although post-depositional dissolution cannot be excluded.

Figure 11Global changes in deep-sea carbonate burial over the last glacial cycle. Curves show, from top to bottom, a global stack of benthic foraminifera δ¹⁸O (Lisiecki and Raymo, 2005); reconstructed bulk sediment MAR; CaCO₃ MAR assuming constant sedimentation rate; CaCO₃ MAR using reconstructed sedimentation rate. The latter three curves show the mean, ±1σ, and ±2σ; see details in the Supplement). Vertical bars highlight cold periods (MIS2/LGM, MIS4, and MIS6) in blue and full interglacial conditions (Holocene and MIS5e) in yellow.

Download

5.2 Changes in deep-sea [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$]

The fact that deep-sea CaCO₃ burial did not increase during the glacial, despite the loss of marginal carbonate depocenters on shelves due to sea level fall, may seem counterintuitive. If the alkalinity input flux from weathering had remained constant (Gibbs and Kump, 1994), the loss of the coastal alkalinity sink should have led to a buildup of ocean alkalinity and a compensatory increase in deep-sea CaCO₃ burial (Opdyke and Walker, 1992). However, this chain of expectations does not consider possible increases in the DIC of the deep sea that could have counterbalanced the impact of an alkalinity increase on the [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] there. Because the [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] is dependent not only on ALK but also on DIC (since [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] ≈ ALK – DIC), large changes in DIC will cause [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] to change in addition to changes caused by variations in ALK. In fact, it seems likely that deep-sea DIC changes would have arisen both from enhanced carbon sequestration by saturation, soft tissue and disequilibrium carbon storage (Eggleston and Galbraith, 2018; Ödalen et al., 2018), and changes in the active carbon inventory (Yu et al., 2013).

As mentioned above, our mass-balance calculations resolve only changes in the total inventories due to geological inputs and outputs, and do not include redistributions of carbon. As a result, they infer changes in [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] that are unrealistically large and disagree with the relatively small changes in [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] that have been reconstructed from foraminiferal proxies (Yu et al., 2013, 2014). We suggest that the main reason for this discrepancy is that large increases of soft tissue and disequilibrium DIC occurred in the deep sea during glacial periods. In support of this, simulations with complex Earth system models driven with glacial boundary conditions have shown increases of both soft tissue and disequilibrium DIC of the same order as the [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] excursions calculated here (order 40 µM each; Eggleston and Galbraith, 2018), which would have been further supplemented by the greater solubility of CO₂ in cold water. Thus, the discrepancy between inferred and observed [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] changes across the last glacial cycle can likely be resolved with the long-suspected changes in the glacial storage of carbon within the deep ocean.

5.3 Important uncertainties in shelf burial

As shown in Fig. 1, modern estimates suggest that approximately half of the global CaCO₃ burial, and most of the C_org burial, occurs on continental margins. Unfortunately, we do not have direct records of changes in shelf burial over time, so we have instead evaluated the likely impact of global shelf burial changes by considering a simple but illustrative scenario, assuming linear relation between shelf surface variation and burial (Fig. 9).

However, organic and inorganic carbon burials have opposing effects on [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$], with the buildup of DIC during low sea level counteracting the buildup of alkalinity, so that the net effect depends on the CaCO₃ : C_org ratio of net burial. For a CaCO₃ : C_org ratio of approximately 1.3 mol mol⁻¹, the opposing effects of DIC and alkalinity perfectly compensate, leading to muted changes of [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] and atmospheric CO₂, even with large burial and inventory changes. This sensitivity to the shelf burial ratio is illustrated in Fig. 12. Given that estimates of the shelf burial ratio range widely and may have changed significantly over time, it does not appear possible to definitively conclude whether the net effect was to raise or lower CO₂ during glacial periods without improved observational constraints.

Figure 12Atmospheric PCO₂, DIC and alkalinity in the deep ocean for the shelf scenarios (Fig. 8b). Each experiment appears as a straight line of dots and each dot corresponds to a 1 kyr time step. The tilt of the individual scenario lines depends on the CaCO₃ : C_org burial ratio in shelf sediments. Blue and red lines correspond to 100 % carbonate and 100 % C_org burial, respectively.

Download

5.4 Implications of geological fluxes for whole-ocean δ¹³C and terrestrial biosphere

Because the changes in burial on shelves would have changed the global CaCO₃ : C_org burial, they must have left a signature in the δ¹³C of the active carbon inventory, with a strong likelihood for lower δ¹³C during MIS3–MIS2. This isotopic effect arises from the fact that C_org burial (which removes carbon with very low δ¹³C) is strongly focused on continental margins, compared to carbonate burial, which is more evenly distributed between margins and the deep sea. Higher C_org burial in the deep sea during glacial maxima (Cartapanis et al., 2016) was insufficient to compensate for reduced burial on shelves. The smaller proportion of low δ¹³C_org buried during periods of low sea level would have produced a downward drift of the δ¹³C of the active carbon inventory (Wallmann et al., 2016), while enhanced volcanic activity could have reduced δ¹³C during deglaciation (Huybers and Langmuir, 2009; Roth and Joos, 2012).

This result, previously identified by Wallmann (2016), is important given the frequently applied technique of calculating changes in the terrestrial biosphere using reconstructed changes in whole-ocean δ¹³C and assuming a constant inventory (Shackleton, 1977; Bird et al., 1996; Ciais et al., 2012; Menviel et al., 2012; Peterson et al., 2014; Schmittner and Somes, 2016). The magnitude of ocean δ¹³C changes that result from changes in shelf burial is highly uncertain, given the large range in estimated modern burial fluxes, but could potentially exceed the observed glacial–interglacial change of 0.32 ‰ (Wallmann et al., 2016). Unfortunately, this invalidates the use of δ¹³C as a direct monitor of changes in the terrestrial biosphere, and requires that new types of constraints on changes in terrestrial biosphere be developed.

Although the uncertainties are large, our mass-balance calculations suggest that the terrestrial biosphere during the LGM was most likely larger than implied by the constant-inventory assumption. We note that using the highest modern C_org burial flux (on the order of 500 PgC kyr⁻¹) would imply a change of the δ¹³C of the system on the order of 1 ‰ (Figs. 8 and 10), which is higher than reconstructed variability for δ¹³C in the deep ocean, in which most of the carbon resides (Oliver et al., 2010; Peterson et al., 2014). This suggests that either the highest estimates for modern burial over shelves (above 500 PgC kyr⁻¹) are inconsistent with the assumption of a 60 % decrease of the shelf burial during peak glacial or that the terrestrial biosphere actually contained more biomass during the glacial than today, in contrast to the standard assumption. Although this is highly uncertain, such a change could have occurred through an expansion of peatlands and permafrost carbon during the glacial (Ciais et al., 2012). An expansion of terrestrial carbon storage in a colder state would perhaps be more consistent with the general expectation for carbon release from peatland and permafrost under anthropogenic warming (IPCC, 2014), though it should be viewed as a highly speculative possibility.

5.5 Non-steady-state Holocene carbon cycle

Our results show that the removal rate of carbon from the active pool was higher during the Holocene (prior to anthropogenic carbon release) compared to the glacial (Fig. 10c). This result, which emerges primarily from the fact that CaCO₃ burial, both on the shelves and in the deep sea, was faster than it had been during the glacial, indicates that the active carbon inventory was shrinking prior to anthropogenic inputs. The estimates suggest that this imbalance was mostly attributable to the shelf burial, though at a rate that was very slow compared to subsequent fossil fuel emissions (∼100 PgC kyr⁻¹ net loss vs. present-day emissions of ∼10 000 PgC kyr⁻¹).

In support of this large glacial–interglacial variation in carbon burial rates, we find that only mass-balance calculations that include reduced burial during the glacial infer background inputs of carbon that are consistent with the <400 PgC kyr⁻¹ geological inputs estimated in the literature (Fig. 1). Essentially, without slower burial during the glacial, modern-day burial estimates would require larger carbon input fluxes than supported by observational estimates in order to balance the long-term inventory. Although this consistency does not amount to proof, it provides encouraging support for the emergent picture of a geological carbon cycle that is dynamic on the timescale of glacial cycles.