Low terrestrial carbon storage at the Last Glacial Maximum: constraints from multi-proxy data

Past changes in the inventory of carbon stored in vegetation and soils remain uncertain. Earlier studies inferred the increase in the land carbon inventory (Δland) between the Last Glacial Maximum (LGM) and the preindustrial period (PI) based on marine and atmospheric stable carbon isotope reconstructions, with recent estimates yielding 300–400 GtC. Surprisingly, however, earlier studies considered a mass balance for the ocean–atmosphere–land biosphere system only. Notably, these studies neglect carbon exchange with marine sediments, weathering–burial flux imbalances, and the influence of the transient deglacial reorganization on the isotopic budgets. We show this simplification to significantly reduce Δland in simulations using the Bern3D Earth System Model of Intermediate Complexity v.2.0s. We constrain Δland to ∼850 GtC (median estimate; 450 to 1250 GtC ±1SD) by using reconstructed changes in atmospheric δ¹³C, marine δ¹³C, deep Pacific carbonate ion concentration, and atmospheric CO₂ as observational targets in a Monte Carlo ensemble with half a million members. It is highly unlikely that the land carbon inventory was larger at LGM than PI. Sensitivities of the target variables to changes in individual deglacial carbon cycle processes are established from transient factorial simulations with the Bern3D model. These are used in the Monte Carlo ensemble and provide forcing–response relationships for future model–model and model–data comparisons. Our study demonstrates the importance of ocean–sediment interactions and burial as well as weathering fluxes involving marine organic matter to explain deglacial change and suggests a major upward revision of earlier isotope-based estimates of Δland.

1 Introduction

Atmospheric CO₂ varied between about 180 and 300 ppm over the past 800 000 years (Neftel et al., 1982; Siegenthaler et al., 2005; Lüthi et al., 2008; Marcott et al., 2014; Bereiter et al., 2015). These CO₂ variations were tightly coupled to glacial–interglacial climate change (Petit et al., 1999; Siegenthaler et al., 2005) and amplified orbitally driven climate variations (Jansen et al., 2007). Despite their importance, the mechanisms behind these past CO₂ variations remain enigmatic. A wide range of explanatory processes related to the marine and terrestrial carbon cycle as well as to exchange processes with reactive ocean sediments, coral reefs, and the lithosphere has been proposed (e.g., Archer et al., 2000; Sigman and Boyle, 2000; Sigman et al., 2010; Fischer et al., 2010; Menviel et al., 2012; Ushie and Matsumoto, 2012; Jaccard et al., 2013; Galbraith and Jaccard, 2015; Heinze et al., 2016; Wallmann et al., 2015). While the first simulations covering glacial–interglacial cycles, including dynamic ocean and land models, represent reconstructed atmospheric CO₂ variability emerge, these models are not yet able to reproduce variations in important proxy data or to represent the timing of CO₂ changes over the last glacial termination adequately (Brovkin et al., 2012; Menviel et al., 2012; Ganopolski and Brovkin, 2017). To make progress in this research area requires the combination of multiple lines of evidence whereby proxy reconstructions are compared to quantitative model analyses. This will help to constrain underlying processes and to quantify their contribution.

One of the processes is the storage of carbon in the land biosphere, but its change over the last deglaciation is debated. Most studies addressing past land biosphere carbon focus on the change in land biosphere carbon inventory between the Last Glacial Maximum around 21 000 years before present (LGM; 21 ka) and the recent preindustrial period (PI). This inventory difference (PI minus LGM) is here termed Δland and includes changes in the various land biosphere reservoirs, such as plants, mineral soils, permafrost, peatland, yedoma, and wetlands as well as changes on shelves exposed during the LGM.

Available proxy-based Δland estimates (e.g., Shackleton, 1977; Adams et al., 1990; Bird et al., 1994; Adams and Faure, 1998; Zimov et al., 2006, 2009; Zech et al., 2011; Ciais et al., 2012; Peterson et al., 2014; Menviel et al., 2017) encompass values ranging from −400 to 1500 GtC. This is large compared to the current land biosphere stock amounting to around 4000 GtC (Ciais et al., 2013) as well as to the deglacial atmospheric change of around 200 GtC (1 ppm is equivalent to 2.12 GtC).

The classical (Shackleton, 1977) and arguably most reliable approach to reconstruct Δland considers the deglacial change in the marine stable carbon isotope signature as recorded in sediments and the mass balance of carbon and ¹³C. The argument is that the uptake of isotopically light carbon by the land biosphere causes a corresponding measurable perturbation in the average δ¹³C value in the ocean–atmosphere system. The ¹³C mass balance approach, using marine sediment and ice core records, provided consistent and converging estimates of Δland with recent estimates ranging from 300 to 400 GtC (Curry et al., 1988; Duplessy et al., 1984; Crowley, 1991; Bird et al., 1994; Crowley, 1995; Bird et al., 1996; Joos et al., 2004; Ciais et al., 2012; Peterson et al., 2014; Menviel et al., 2017). Early estimates relied on spatially limited δ¹³C sediment records (Duplessy et al., 1984; Curry et al., 1988; Bird et al., 1994, 1996), but recently, more comprehensive δ¹³C data compilations (Oliver et al., 2010; Peterson et al., 2014) have become available. These were also used to determine mean ocean δ¹³C and Δland by applying a dynamic ocean model (Menviel et al., 2017).

All previous ¹³C mass balance studies assumed a closed system. They neglected the carbon and isotopic exchange between the ocean and reactive marine sediments, the input flux from rock weathering, and carbon burial in consolidated sediments by focusing on the atmosphere–ocean–land system only. It has been argued (e.g., Ciais et al., 2012) that changes in weathering contributed little to the deglacial CO₂ rise and that large changes in the calcium carbonate (CaCO₃) burial in sediments (and coral reefs) are inconsistent with the reconstructed changes in the depth of the lysocline.

There are, however, several important reasons to question these arguments and the closed system assumption. First, a continuous flux of calcareous biogenic particles and organic matter carries any perturbation in δ¹³C in the coupled atmosphere–surface ocean system to ocean sediments (Tschumi et al., 2011; Menviel et al., 2012; Roth et al., 2014; Cartapanis et al., 2016). A second reason relates to the chemical buffering of land-induced CO₂ perturbations to marine CaCO₃ sediments, a process known as carbonate compensation (e.g., Archer et al., 2000; Broecker and Clark, 2001). Carbonate compensation causes a net transfer of carbon and ¹³C between the ocean and reactive marine sediments. Third, and of even greater importance with regard to δ¹³C, are changes in the cycling and weathering–sedimentation imbalances of organic carbon that have the potential to exert a large impact on the ¹³C budget (Tschumi et al., 2011; Roth et al., 2014). These mechanisms should not be neglected when addressing the whole-ocean budgets of carbon and δ¹³C to infer changes in land carbon stocks.

The transfer of ¹³C and carbon between the ocean–atmosphere system and reactive sediments and the lithosphere is affected by the well-documented reorganization of the marine carbon cycle across the deglaciation (e.g., Sigman and Boyle, 2000; Sigman et al., 2010; Fischer et al., 2010; Elderfield et al., 2012; Ciais et al., 2013; Martinez-Garcia et al., 2014; Jaccard et al., 2016; Cartapanis et al., 2016). The reorganization includes changes in the following: CO₂ solubility, ocean ventilation (e.g., Burke and Robinson, 2012; Sarnthein et al., 2013; Skinner et al., 2017), and water mass distribution (e.g., Duplessy et al., 1988; Peterson et al., 2014; Lippold et al., 2016; Menviel et al., 2017; Gottschalk et al., 2018); air–sea gas transfer rates, for example via changes in sea-ice extent (e.g., Stephens and Keeling, 2000; Gersonde et al., 2005; Waelbroeck et al., 2009; Sun and Matsumoto, 2010), export (e.g., Kohfeld et al., 2005; Jaccard et al., 2013), and the remineralization (e.g., Bendtsen et al., 2002; Matsumoto, 2007; Taucher et al., 2014) of biogenic material associated with the cycling of organic carbon, CaCO₃, and biogenic opal (Tschumi et al., 2011; Roth et al., 2014; Matsumoto et al., 2014); coral reef growth (e.g., Berger, 1982; Milliman and Droxler, 1996; Kleypas, 1997; Ridgwell et al., 2003; Vecsei and Berger, 2004; Menviel and Joos, 2012) and the input of dust, nutrients, and lithogenic material by atmospheric deposition (e.g., Lambert et al., 2015) and rivers; and coastal and shelf exchange processes (Ushie and Matsumoto, 2012; Wallmann et al., 2015). These changes influence the particle flux of organic and inorganic carbon, opal, and mineral particles, adding carbon and mass to the sediments. As a consequence, these deglacial changes and their spatiotemporal evolution affect the mean δ¹³C signature of the global ocean and estimates of Δland. For this reason, it is essential to explicitly consider the causes of glacial–interglacial variations in CO₂ to quantify Δland.

The primary goal of this study is to provide a new, observationally constrained best estimate of Δland recognizing the role of marine deglacial processes, weathering and burial, and sediment interactions in an Earth system model. To this end, we establish the response sensitivities for mean ocean δ¹³C_DIC and for other proxy targets (PI–LGM differences) to the changes in individual key deglacial carbon cycle mechanisms from transient simulations. An emulator of the Bern3D has been developed, and the strengths of individual marine and terrestrial processes are varied in a large number of combinations to find all possible solutions for Δland that match a set of observational targets in a Bayesian Monte Carlo data assimilation framework (see Fig. 1). We further scrutinize the closed system assumption in the ¹³C mass balance approach with land carbon uptake as the only forcing. Specifically, we run idealized pulse-like carbon uptake simulations with the Bern3D Earth System Model of Intermediate Complexity (EMIC) both with and without simulating sediment interactions and with and without interactive weathering. We further discuss multi-proxy response relationships for different carbon cycle processes and briefly address the possible contribution of individual mechanisms to the deglacial change in atmospheric CO₂ and δ¹³C_atm and in the deep ocean carbonate ion concentration. Our results demonstrate the importance of different deglacial mechanisms, ocean–sediment interactions, and the burial–weathering cycle for the whole-ocean δ¹³C_DIC signature. All previous δ¹³C-based estimates of Δland neglected the influence of these processes and appear systematically biased towards lower estimates.

Figure 1 Flow chart outlining the steps applied to achieve an estimate of possible Δland values. To this end, seven generic deglacial carbon cycle mechanisms (Table 2 and Sect. 2.2) and four observational targets (Sect. 2.2.3) were included. The seven processes are varied individually by systematic parameter variations in addition to the well-established forcings (top of figure) in 50 factorial simulations with the Bern3D. This yields a first set of sensitivities or forcing–response curves for our four targets and seven processes. Next, we apply a Latin hypercube parameter sampling to vary the processes in combination and probe for nonlinear interactions in 60 multiparameter simulations with the Bern3D model. The results are used to adjust the forcing–response relationship; adjustments are small, except for CO₂ response curves (see Fig. C1 in the Appendix). Analytical representations of the forcing–response curves are used to build a simple emulator of the Bern3D model (Sect. 2.2.4). The emulator is applied in a Monte Carlo ensemble with half a million members to probe a wide range for each individual process. Finally, the ensemble results are constrained by the four proxy targets to yield a best estimate and an uncertainty distribution for Δland.

2 Methods

2.1 The Bern3D model

The Bern3D EMIC couples a dynamic geostrophic-frictional balance ocean, a thermodynamic sea-ice component, and a single-layer energy–moisture balance atmosphere. The horizontal resolution is 41×40 grid cells and the ocean has 32 layers. Marine productivity is simulated as a function of nutrient concentrations (P, Fe, Si), temperature, and light and transferred to dissolved organic and particulate matter. Biogenic matter decomposes within the water column. Opal, calcite, and organic particles reaching the ocean floor are entering reactive sediments. A 10-layer sediment model is used to compute fluxes of carbon, nutrients, alkalinity, and isotopes between the ocean, reactive sediments, and the lithosphere. Loss fluxes to the lithosphere are compensated for at equilibrium by a corresponding input flux from weathering. A four-box reservoir model is used to calculate the dilution of atmospheric isotopic perturbations by exchange with the land biosphere.

Further information on the model and the model spin-up and experimental settings for the impulse response experiments is given in the Appendix A. The sediment module has been updated to include the most recent observational information on sediment composition and fluxes (Cartapanis et al., 2016, 2018; Tréguer and De La Rocha, 2013) as detailed in Appendix B. Key model parameters were adjusted within their uncertainties to best reproduce the observation-based spatial distributions and global stocks of CaCO₃, particulate organic carbon (POC), and opal within reactive sediments, observational estimates of global burial and redissolution fluxes, and euphotic-zone export fluxes of biogenic particles (Table 1). This updated version of the Bern3D model is introduced as Bern3D v2.0s. The simulated oceanic dissolved inorganic carbon (DIC) inventory for PI is 37 175 GtC, similar to the 37 310 GtC estimated based on the GLODAP v.2 dataset (Lauvset et al., 2016).

Figure 2 Observed (a) versus simulated (b) preindustrial δ¹³C_DIC distributions along a cross section through the Atlantic (25^∘ W), across the Southern Ocean (58^∘ S), and into the Pacific (175^∘ W). (a) Data from Eide et al. (2017) based on ocean measurements and corrected for the Suess effect. (b) δ¹³C as modeled in the Bern3D model under preindustrial boundary conditions.

¹³C stocks and exchanges are modeled in the atmosphere–ocean–sediment–land biosphere system. ¹³C fluxes to the lithosphere associated with the burial of POC and CaCO₃ and weathering input to the ocean are explicitly simulated. The burial of POC and CaCO₃ results in an isotopic signature of the burial flux of around −9 ‰, intermediate between the isotopically light POC ($\sim -\mathrm{20}$ ‰) and the isotopically heavier CaCO₃ signature (∼2.9 ‰; see Appendix A for implementation of isotopic discrimination). An isotopic perturbation in the atmosphere, e.g., through uptake of isotopically light carbon by the land biosphere, is transferred to the ocean through air–sea gas exchange and to the land by photosynthesis. Within the ocean, any perturbation is communicated between the surface ocean and the deep ocean by advection, convection, and mixing. In addition, POC and CaCO₃ export fluxes from the surface ocean contribute to the burial flux.

Figure 2 shows modeled δ¹³C in a section through the Atlantic, Southern Ocean, and Pacific under preindustrial (1765 CE) boundary conditions with a prescribed atmospheric δ¹³C signature of −6.305 ‰. Compared to measured data corrected for the Suess effect (Eide et al., 2017) it is visible that modeled δ¹³C values are biased towards high values by a seemingly constant offset of about 0.4 ‰, while the δ¹³C patterns are faithfully captured by the model.

Table 1 Export, burial, and deposition fluxes, as well as sediment inventories of CaCO₃, opal, and particulate organic carbon (POC) as determined in the model after a preindustrial spin-up and literature estimates. The tracer input flux to the ocean resulting from weathering is set to compensate for the burial fluxes diagnosed at the end of the spin-up and kept constant in Bern3D standard simulations.

^a Battaglia et al. (2016), ^b Tréguer and De La Rocha (2013), ^c Sarmiento and Gruber (2006), ^d Milliman and Droxler (1996), ^e Feely et al. (2004).

2.2 Transient LGM to PI sensitivity simulations

2.2.1 Standard transient forcings

Starting from PI steady-state conditions, the model is first run for 20 kyr under constant LGM forcing and then for an additional 20 kyr under transient forcing from the LGM to PI. Forcings (Fig. 3) include radiative forcing imposed by CO₂, N₂O, and CH₄ (Joos and Spahni, 2008) and variations in orbital parameters (Berger, 1978). Radiative forcing of CO₂ is prescribed and biogeochemistry is determined with interactive CO₂. Ice sheet extent and related changes in albedo are prescribed based on benthic δ¹⁸O data (Lisiecki and Raymo, 2005) and ice sheet reconstructions (Peltier, 1994). Additionally, freshwater (FW) pulses are prescribed in the North Atlantic as detailed in Menviel et al. (2012), also accounting for the LGM to PI change in salinity of about 1 PSU. Coral reef regrowth is implemented from 14 ka onward (Vecsei and Berger, 2004) by uniformly removing the respective amount of alkalinity, carbon, and ¹³C from the uppermost ocean grid cells. These forcings (Fig. 3a–d) are termed standard forcings.

Figure 3 Forcings applied for the LGM to PI simulations. (a) Radiative forcing imposed by CO₂, CH₄, and N₂O, (b) benthic δ¹⁸O used to scale ice sheet extent, (c) freshwater forcing into the North Atlantic, and (d) coral reef regrowth as applied as standard transient forcing. (e) Prescribed evolution of land biosphere carbon uptake. Cumulative land biosphere uptake is varied over the termination across factorial simulations but kept invariant (254 GtC) over the Holocene (11–0 ka) following Elsig et al. (2009); here a scenario with a total uptake of 445 GtC is shown for illustration. (f) Idealized evolution of scaling factors as prescribed in factorial experiments to vary the remineralization profile, rain ratio, Southern Ocean wind stress and gas transfer rate, or organic weathering flux on land between LGM and PI values. The freshwater forcing in the North Atlantic is varied over 40 kyr. Dashed green lines indicate spin-up values.

2.2.2 Factorial sensitivities

Standard transient deglacial forcings in the Bern3D model are complemented by changes in seven additional mechanisms, grouped into generic classes and introduced below (see top of Fig. 1). A set of 50 factorial experiments (Table 2), covering the same time interval as the model run with standard transient forcings, is conducted to quantify the sensitivities of the carbon cycle to changes in these seven mechanisms. The forcing history for changes in Southern Ocean wind stress, Southern Ocean gas transfer rate, the rain ratio, the remineralization profile, and the organic weathering rate are prescribed following an idealized glacial–interglacial evolution. PI values are set to LGM values at 40 ka, kept constant from 40 to 18 ka, scaled back to the PI value over the termination (18 to 11 ka), and kept constant thereafter (Fig. 3f and Table 2). The sensitivities of carbon cycle properties are obtained by subtracting the standard forcing run from the runs with complementary changes in the seven mechanisms. This also yields characteristic response relationships between different carbon cycle properties, e.g., changes in atmospheric CO₂ versus changes in whole-ocean δ¹³C, for each process. The objective is to schematically represent the space of plausible property–property relationships for key proxies. In other words, the selected mechanisms are assumed to be representative on the global scale for the myriad of processes influencing the carbon cycle and will be outlined below.

Physical mechanisms. Changes in ocean circulation, for example in response to altered wind stress, stratification, eddy mixing, or buoyancy forcing, affect the marine carbon inventory (e.g., Adkins, 2013) by altering the cycling of organic and inorganic carbon. The Southern Ocean has been identified as an important region modulating CO₂ changes, and a reduced Antarctic overturning circulation has been invoked to, at least partly, explain low glacial atmospheric CO₂ (Sigman and Boyle, 2000; Sigman et al., 2010; Fischer et al., 2010; Skinner et al., 2010; Tschumi et al., 2011; Ferrari et al., 2014; Watson et al., 2015; Roberts et al., 2016; Gottschalk et al., 2016; Jaccard et al., 2016). Under standard forcing, Southern Ocean circulation changes are weak in the Bern3D model. In this study, the circulation is modified by uniformly scaling the wind stress field over the Southern Ocean ($>{\mathrm{48}}^{\circ }$ S, south of South America) (see also Tschumi et al., 2008) to reduce deep ocean ventilation. We do not consider such large changes in wind stress to be realistic, but rather use prescribed Southern Ocean wind stress as a way to vary deep ocean ventilation in the Bern3D model. Consequently, we do not scale air–sea gas transfer rates with the changes in wind stress but vary transfer rates independently.

Proxy reconstructions suggest a larger sea-ice extent in the Southern Ocean during glacial times (Gersonde et al., 2005; Waelbroeck et al., 2009). Sea ice hinders air–sea gas exchange and thus exerts a large impact on atmospheric δ¹³CO₂ (Lynch-Stieglitz et al., 1995). Additionally, changes in Southern Ocean winds (Kohfeld et al., 2013) may have affected glacial Southern Ocean air–sea gas transfer rates. Here, the standard air–sea gas transfer rate is scaled (Fig. 3f; Table 2) in the Southern Ocean to estimate related carbon cycle sensitivities. We note that changes in sea ice and a corresponding change in air–sea gas transfer rates are reasonably well simulated by the Bern3D model. The effect of changes in ocean temperature (Waelbroeck et al., 2009) and salinity on the solubility of CO₂ is discussed elsewhere (Sigman and Boyle, 2000; Brovkin et al., 2007; Sigman et al., 2010; Menviel et al., 2012) and is part of the standard forcing.

Table 2 Overview of mechanisms and parameters considered in factorial sensitivity experiments. Bold entries mark the standard parameter values. Values represent the export rain ratio at LGM, the change in remineralization profile (see text) at LGM, the scaling applied to the standard shallow water carbonate deposition history, and the land carbon uptake over the deglacial in gigatons of carbon (GtC). In the case of SO wind stress, SO gas transfer rate, and the organic weathering flux, values represent the ratio of LGM to PI forcing.

Oceanic carbonate. The partial pressure of CO₂ in the ocean depends on alkalinity and processes that alter surface ocean alkalinity, thus influencing atmospheric CO₂. The sedimentary burial of CaCO₃ removes twice as much alkalinity as carbon from the ocean, leading to oceanic CO₂ outgassing. Increased CaCO₃ burial over the deglacial has been put forward as the coral reef hypothesis (Berger, 1982) to explain increases in atmospheric CO₂. Vecsei and Berger (2004) reconstructed coral reef growth history and provided a lower-limit estimate of about 380 GtC for the amount of shallow water carbonate deposition over the deglacial period. Other estimates of CaCO₃ deposition amount to 1200 GtC or even more (Milliman, 1993; Kleypas, 1997; Ridgwell et al., 2003). This yields a large range of possible scenarios with substantial impacts on atmospheric CO₂. Here, we scale the reconstructed deposition history based on Vecsei and Berger (2004) (Fig. 3d) by applying a constant scaling to vary CaCO₃ deposition within the published range.

Potential changes in the rain ratio (CaCO₃ : C_org) affect both alkalinity and CO₂ (e.g., Berger and Keir, 2013; Archer and Maier-Reimer, 1994; Sigman et al., 1998; Archer et al., 2000; Matsumoto et al., 2002; Tschumi et al., 2011). An increase in the export of CaCO₃ relative to C_org from surface waters has similar consequences for CO₂ as shallow water carbonate deposition (Sigman and Boyle, 2000). There is a lack of estimates on how the rain ratio has varied in the past. Here, we vary the global rain ratio for LGM conditions (see also Appendix A; Fig. 3f; Table 2) between 0.045 and 0.098, while the rain ratio is 0.083 in the standard setup and for PI conditions.

Oceanic organic matter. There is a broad range of mechanisms, in addition to circulation changes, that affect the cycling of organic matter within the ocean, the whole-ocean nutrient inventories, and thereby surface ocean nutrient and DIC concentrations as well as atmospheric CO₂. Several studies have suggested changes in the remineralization length scale of organic carbon due to colder ocean temperatures (see, e.g., Matsumoto, 2007; Menviel et al., 2012; Roth et al., 2014). A deeper remineralization of organic matter leads to the removal of nutrients and DIC from the shallow subsurface ocean and, via ocean–sediment interactions, to globally reduced nutrient and increased alkalinity inventories (Menviel et al., 2012; Roth et al., 2014), resulting in a decrease in atmospheric CO₂ (Matsumoto, 2007; Menviel et al., 2012). Weathering–burial feedbacks lead to an amplification of changes in atmospheric and oceanic δ¹³C as shown by Tschumi et al. (2011) and Roth et al. (2014). Here, the remineralization profile in the upper 2 km is varied between the standard Martin curve (Eq. A1 in Appendix) and a linear profile (75 to 2000 m). The modification in the remineralization profile of particulate organic matter is different than those applied by Menviel et al. (2012) and Roth et al. (2014) in simulations with the Bern3D model. Here, only the remineralization profile in the upper 2000 m is changed, and, in contrast to these earlier studies, it is assumed that the fraction of exported particles that reach the sediments below 2000 m is not altered. The implicit assumption is that cooler temperatures at the LGM only significantly affect the dissolution of relatively small or labile particles in the upper ocean, while large or refractory particles sink fast enough to reach ocean sediments both under LGM and PI conditions.

Oceanic organic carbon and nutrient inventories are further affected by changes in the weathering flux compensating for the burial of organic matter (Cartapanis et al., 2016, 2018). A negative balance between the riverine input of PO₄, DIC, DI¹³C, and Alk, originally buried as organic material, and sedimentary burial of particulate organic matter leads to a transient decrease in atmospheric CO₂ and also affects atmospheric and marine δ¹³C. The weathering flux to the ocean (Table 1) is kept constant in the standard setup of the Bern3D (see also Appendix A) but might have varied with changes in climate and CO₂. To account for this possibility, we varied the input of PO₄, DIC, DI¹³C, and Alk that compensates for the burial of particulate organic material in the sediment in factorial experiments (Fig. 3f; Table 2). This variation in the input will be referred to as changes in the organic weathering flux.

Another widely discussed mechanism relates to iron fertilization in the Southern Ocean (Martin, 1990; Archer et al., 2000; Sigman and Boyle, 2000; Parekh et al., 2008; Sigman et al., 2010; Fischer et al., 2010; Martinez-Garcia et al., 2014). Iron fertilization leads to changes in the cycling of marine organic material and to related changes in the whole-ocean inventory of δ¹³C_DIC and other variables. These changes are broadly similar to the changes caused by altering the organic matter remineralization depth or the weathering input flux and implicitly included in the Monte Carlo simulations. Similarly, other processes, such as changes in the ratio of nutrient to carbon uptake by marine organisms (Heinze et al., 2016) or changes in the association of POC with ballast minerals (Armstrong et al., 2002), may have affected the efficiency of the marine biological cycle in reducing carbon in the surface ocean and the atmosphere.

Land biosphere carbon inventory. Changes in land biosphere carbon stock are prescribed following the evolution illustrated in Fig. 3e. We assume that the change in land carbon stock over the Holocene is well constrained and amounts to an uptake of 290 GtC at 11–5 ka and a release of 36 GtC thereafter (Elsig et al., 2009). Scenarios with different cumulative land carbon change imply a corresponding uptake or release of terrestrial carbon during the last glacial termination (18–11 ka). The δ¹³C signature of terrestrial carbon is set to −24 ‰.

2.2.3 Data constraints for deglacial carbon and carbon isotope changes

We consider four globally relevant observational constraints (Table 3) to estimate the PI–LGM change in the land biosphere carbon inventory, Δland. These are the PI–LGM difference (Δ) in (i) atmospheric CO₂ (ΔCO₂), (ii) atmospheric δ¹³C(CO₂) (Δδ¹³C_atm), (iii) mean ocean δ¹³C(DIC) (Δδ¹³C_DIC), and (iv) deep (equatorial) Pacific carbonate ion concentration ($\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$) (Table 3). The carbon isotope constraints are directly relevant for the isotopic mass balance used to infer Δland. The carbonate ion concentration difference constrains changes in the ocean's alkalinity budget, and the change in atmospheric CO₂ constrains the integrated effect of deglacial carbon cycle changes. We allow for a relatively wide uncertainty range of 20 ppm for ΔCO₂, which is admittedly wider than the proxy uncertainty, to avoid an overfitting of model outputs. The Δδ¹³C_atm range for the constraint is 0.1±0.05 ‰ and based on ice core measurements (Schmitt et al., 2012). Peterson et al. (2014) compiled 480 benthic foraminiferal records and report a whole-ocean PI minus LGM change in δ¹³C_DIC of 0.34±0.19 ‰, and we use their range (0.15 ‰ to 0.53 ‰) as a constraint. Proxy reconstructions of $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ based on measured B∕Ca ratios in foraminifera show little change from LGM to PI in the equatorial deep Pacific (Yu et al., 2013). Based on a different approach, Qin et al. (2017) also conclude on negligible change in $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ in the western tropical Pacific across the last glacial termination. In a recent study, Luo et al. (2018) report, based on CaCO₃ reconstructions in the South China Sea, that carbonate ion concentrations in the mid-depth Pacific hardly changed from LGM to PI. All these studies point to a small change in $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$, providing a further constraint on deglacial simulations. However, the uncertainty in these reconstructions is of the same order as the change itself. The target range for $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ is chosen here from −6 to 5 µmol kg⁻¹.

Schmitt et al. (2012)Peterson et al. (2014)

Table 3 Overview of the four proxy-based constraints representing change (Δ) between PI minus LGM. LGM refers to the period 20 to 19 ka for the ice core and to the period 23 to 19 ka for the ocean sediment data. Similarly, PI refers to 500 to 200 BP for the ice core data and 6000 to 200 BP and 8 to 6 ka for marine δ¹³C and ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$, respectively. The assumed uncertainty range in ΔCO₂ is considered to be wider than the uncertainty from the ice core data. This is to avoid overfitting and taking into account uncertainties in the emulator introduced in Sect. 2.2.4.

2.2.4 Emulator and Monte Carlo setup

Sensitivities of the four target proxies (ΔCO₂, Δδ¹³C_atm, Δδ¹³C_DIC, and $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$) to changes in the selected seven forcings are computed based on 50 factorial experiments. For each forcing, i, and related parameter change, Δp_i (Table 2), the corresponding sensitivity, i.e., the change in target per unit change in parameter, is computed. These discrete sensitivities are fitted linearly or quadratically for each target, T, and forcing to get a continuous function for the sensitivity, $S_i^T$(Δp_i). These sensitivity functions are used to build a simple representation to compute changes in target variables for a combination of changes in parameters:

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }T\left(\mathrm{\Delta }{p}_{\mathrm{1}},\mathrm{\dots },\mathrm{\Delta }{p}_{\mathrm{7}}\right)=\sum _{i=\mathrm{1}}^{\mathrm{7}}\mathrm{\Delta }{p}_i\times S_i^T\left(\mathrm{\Delta }{p}_i\right).\end{array}$$

In order to evaluate the assumed linear additivity of the seven mechanisms, a 60-member ensemble of parameter combinations is calculated both with the Bern3D model and Eq. (1) (Fig. 1). The seven parameters are sampled uniformly using Latin hypercube sampling (LHS; McKay et al., 1979). The results from the Bern3D model and Eq. (1) are regressed against each other for each target variable (Fig. C1). The corresponding linear fit is used to account for nonadditive behavior when combining forcings. This then yields the following emulator:

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }T\left(\mathrm{\Delta }{p}_{\mathrm{1}},\mathrm{\dots },\mathrm{\Delta }{p}_{\mathrm{7}}\right)=a^T+b^T\times \left(\sum _{i=\mathrm{1}}^{\mathrm{7}}\mathrm{\Delta }{p}_i\times S_i^T\left(\mathrm{\Delta }{p}_i\right)\right),\end{array}$$

where a^T is the offset and b^T the slope of the respective linear fit as given in Appendix Fig. C1. The root mean square deviation between the emulator (Eq. 2) and the Bern3D model is 8.8 ppm for ΔCO₂, 0.06 ‰ for Δδ¹³C_atm, 0.02 ‰ for Δδ¹³C_DIC, and 4 µmol kg⁻¹ for $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$. The seven parameters of the LGM to PI sensitivity experiments are sampled uniformly by drawing random samples over their employed ranges (Table 2) to generate a 500 000-member Monte Carlo ensemble. The above emulator is used for each of these members to estimate the corresponding changes in the four targets. Members that yield changes in the targets that lie outside the defined ranges (Table 3) are excluded from the analysis; all remaining members are used to calculate median and confidence ranges in Δland and other variables of interest.

Figure 4 Temporal evolution of the perturbation in (a) atmospheric CO₂ and δ¹³C_atm. (b) The calculated perturbation in terrestrial carbon storage following Eqs. (A2) and (A3) and δ¹³C_DIC anomalies for a pulse uptake of 100 GtC by the land in a closed (green), open (orange), and open system with enabled weathering feedback (purple). Solid lines refer to left and dashed lines to right y axis in both panels.

3 Results

3.1 Pulse experiments

We start by investigating how atmospheric CO₂ and δ¹³C signatures evolve in response to carbon uptake by the land biosphere when all other forcings remain unchanged. To this end, 100 GtC of light carbon (${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}=-\mathrm{24}$ ‰) is instantaneously (in the first time step of the year) removed from the model atmosphere to test the Bern3D response for different model setups (see Appendix A). This yields the Green's functions (or impulse response functions) for carbon and ¹³C shown in Fig. 4.

First, we apply the pulse to a closed system, which includes only the atmosphere, ocean, and land biosphere components. The removal of 100 Gt of light carbon from the atmosphere results in an initial decline in atmospheric CO₂, followed by a typical impulse–response recovery (e.g., Maier-Reimer and Hasselmann, 1987; Siegenthaler and Joos, 1992). Atmospheric CO₂ recovers as the ocean starts to outgass CO₂ (Fig. 4). In a closed system, atmospheric CO₂ reaches a new equilibrium after about 2 kyr, at 8 ppm below the initial steady state. As light carbon is removed, the remaining carbon in the ocean and atmosphere becomes relatively enriched in ¹³C. Initially, the δ¹³C_atm perturbation is removed much faster than the atmospheric CO₂ perturbation as exchange with the ocean and land biosphere dilutes the imposed isotopic perturbation (Fig. 4a; Joos and Bruno, 1996). The oceanic δ¹³C increases and equilibrates after about 2 kyr. The difference in the response of CO₂ and δ¹³C is related to their different properties. CO₂ represents a concentration, whereas δ¹³C stands for the isotopic ratio of ¹³C to ¹²C. The capacity of the ocean to remove an atmospheric perturbation in CO₂, and equally in ¹²CO₂ and ¹³CO₂, is limited by the acid–base carbonate chemistry described by Revelle and Suess (1957). In contrast, the removal of a perturbation in the isotopic ratio ¹³CO₂∕¹²CO₂ is hardly affected by this chemical buffering; the buffering affects ¹³CO₂ in the nominator and ¹²CO₂ in the denominator about equally, leading to a near cancellation of its effect on the ratio. The resulting change for the new equilibrium in δ¹³C amounts to 0.066 ‰ for the atmosphere, 0.060 ‰ for the ocean, and 0.065 ‰ for the four-box land biosphere. The mass balance between ocean, land, and atmosphere is maintained such that the inferred land perturbation (calculated using Eqs. A2 and A3) corresponds to the prescribed removal of 100 GtC (Fig. 4b).

In the second open system experiment, marine sediments are included. The evolution of atmospheric CO₂ (see Archer et al., 1998) is similar to the closed system except that CaCO₃ compensation leads to a further recovery of atmospheric CO₂ until equilibrium is approached with an e-folding timescale of about 14 000 years. CO₂ equilibrates at around 3.5 ppm below the initial steady state. For the first thousand years, the evolution of δ¹³C in the atmosphere and ocean is comparable to the closed system experiment but differs thereafter.

The novel and most important finding from the pulse experiment is that the δ¹³C perturbation in the ocean–atmosphere–land biosphere system is increasingly and eventually completely removed by open system processes. Figure 4b shows the δ¹³C_DIC perturbation to decrease on a multi-millennial timescale in contrast to the closed system assumption (dashed lines in Fig. 4b). In other words, benthic foraminifera record a much smaller perturbation in seawater δ¹³C than expected from closed system modeling. The resulting difference in the δ¹³C_DIC perturbation between the open and closed system causes an erroneous mass balance inference for the terrestrial biosphere when applying conventional closed system equations (Eq. A2 and A3). On a timescale of 10 kyr, underestimation of the terrestrial carbon inventory change amounts to 30 %. The error increases as time evolves.

In detail, the difference between the open and closed system pulse responses is explained as follows. The ¹³C budget in the open system changes due to imbalances between input and burial fluxes of POC and CaCO₃ and due to changes in the δ¹³C signatures of the respective fluxes. The uptake of isotopically light land carbon from the atmosphere leads to an increase in the δ¹³C signature of surface waters and through that in the δ¹³C signatures of exported and eventually buried POC and CaCO₃. The mean δ¹³C signature of the burial flux amounts to −8.7 ‰ (δ¹³C${}_{{\mathrm{CaCO}}_{\mathrm{3}}}$: 2.9 ‰) and C_org (δ¹³C_POC: −20.2 ‰) for the first 10 kyr compared to −9.1 ‰ for the weathering input. This difference tends to mitigate the δ¹³C_DIC perturbation. Further, carbon–climate feedbacks cause changes in temperature and ocean circulation, which affect the export of POC and CaCO₃. Burial fluxes of POC further depend on O₂ concentrations and burial fluxes of CaCO₃ on ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ concentrations that evolve over the course of the simulation. As a result, the cumulative sedimentation–weathering imbalance amounts to a removal of 48.5 GtC during the first 10 kyr. The contribution to differences in the carbon isotopic budget relative to the closed system is dominated by changes in POC cycling and its associated δ¹³C signature, while changes in the CaCO₃ cycle play a smaller but non-negligible role. Overall, ¹³C is lost from the atmosphere–ocean system, and δ¹³C_DIC is decreasing.

In the third pulse experiment, changes in global mean air temperature and atmospheric CO₂ concentration affect weathering fluxes of CaCO₃ and CaSiO₃ following Colbourn et al. (2013). The implemented feedbacks cause a reduction in modeled weathering fluxes in response to decreased atmospheric CO₂ and temperature over the course of the experiment. Decreased weathering fluxes of CaCO₃ and CaSiO₃ lead to a decline in alkalinity supplied to the ocean, resulting in a net transfer of carbon to the atmosphere. This further reduces the initial CO₂ perturbation (Fig. 4a). The mean δ¹³C value related to weathering is −9.1 ‰, thereby leading to a slightly higher δ¹³C_DIC perturbation at 10 kyr compared to the open system experiment with fixed weathering fluxes. The error related to the calculated land perturbation is comparable on a 10 kyr time horizon and increases less as time progresses further compared to the open system setup with fixed weathering.

The pulse responses make it clear that open system processes such as sediment interactions, burial, and weathering affect the carbon and δ¹³C budgets on deglacial timescales. Taken at face value, the results suggest that the classical δ¹³C-based estimates of Δland are substantially biased low because of the neglect of open system processes. However, uptake by land carbon is the only forcing in these idealized pulse simulations, whereas climate and biogeochemical cycles underwent a massive reorganization over the glacial–interglacial transition. This reorganization affected open system processes and the ¹³C budget in addition to carbon uptake by the land biosphere. In the next section, we estimate Δland by taking into account deglacial processes.

Figure 5 Histogram (blue bars, normalized) and kernel density estimate (blue line) of Δland wherein parameter combinations fulfill all four data-based constraints. The thick vertical black line indicates the median Δland of ∼850 GtC, and thin black lines show the ±1σ range. The remaining four lines show kernel density estimates of Δland when considering one data-based constraint at a time.

3.2 Ensemble approach to constrain Δland

We represent the deglacial reorganization of the ocean in the standard model setup, complemented by the seven deglacial carbon cycle mechanisms and constrained by the four observational targets within an open system model framework. Applying the Bern3D emulator (Eq. 2) to our 500 000-member parameter ensemble, we find that the prior range in Δland, ranging from −500 to 1500 GtC, is constrained to ∼450 to ∼1250 GtC (±1σ range; see Fig. 5). The constrained median amounts to ∼850 GtC, and almost no ensemble members that fulfill all four observational constraints are found for negative Δland values. The constrained ensemble thus clearly suggests a positive Δland.

Next, we show that it is not sufficient to use Δδ¹³C_DIC, or any other target, in isolation to constrain Δland. The probability distribution of Δland covers the whole range of sampled values almost uniformly (Fig. 5) when applying only one out of the four targets. Considering either Δδ¹³C_DIC or Δδ¹³C_atm together with ΔCO₂ as constraints moves the Δland distribution towards the multi-proxy-based distribution (blue line in Fig. 5; not shown). We conclude that the incorporation of multiple constraints is essential to narrow uncertainties in Δland.

Figure 6 Histogram (thin lines, normalized) and kernel density estimate (thick lines) of (a) Δδ¹³C_DIC, (b) Δδ¹³C_atm, (c) ΔCO₂, and (d) $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ values for parameter combinations calculated with the emulator that fulfill the observational constraints. Different colors in all four panels correspond to subsamples divided along Δland values that are smaller than 1500, 1000, 750, 450, and 0 GtC. Darker grey shading indicates the range of the respective observational constraint.

Our unconstrained ensemble yields a very large spread in the four target variables. Δδ¹³C_DIC ranges from −0.77 ‰ to 1.25 ‰, Δδ¹³C_atm from −0.9 ‰ to 1.32 ‰, ΔCO₂ from −107 to 181 ppm, and $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ from −49 to 42 µmol kg⁻¹. Within the optimization procedure, we specify that only defined, observation-based ranges of these target variables are allowed (Sect. 2.2.3) and that all four of these constraints have to be met at the same time. This yields specific distributions of the target variables within our constrained ensemble (Fig. 6, black lines). The distribution of Δδ¹³C_atm values, for instance, is almost uniform within the defined observational range (Fig. 6b), while ΔCO₂, $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$, and Δδ¹³C_DIC values cluster towards the lower end of the specified target ranges (Fig. 6c, d, a). This is a result of the underlying carbon cycle–climate processes as tested within Bern3D. For example, it remains difficult to model a PI–LGM ΔCO₂ of ∼90 ppm under the constraint of meeting the other three observational targets. Most ensemble members represent a ΔCO₂ of ∼82 ppm (Fig. 6c). Higher ΔCO₂ values are present in the unconstrained ensemble but are generally associated with even more negative $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ (not shown), which is outside the defined target range for $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$.

We stratify the results by Δland and define four sets in which Δland is allowed to be maximally 1000, 750, 450, and smaller than 0 GtC (Fig. 6, blue, green, orange, and magenta lines). Only one target variable, Δδ¹³C_DIC, is sensitive to this restriction within the constrained ensemble and positive Δland. The distributions of Δδ¹³C_atm, ΔCO₂, and $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ are almost identical for all subsets of Δland with the exception of Δland, which is smaller than 0 GtC. Hence, these three targets do not discriminate per se between smaller and larger estimates of positive Δland. In contrast, restricting Δland to smaller values shifts the realized Δδ¹³C_DIC towards lower values, and the mean of the constrained distribution shifts away from the mean observational estimate of 0.34 ‰ (Fig. 6a). For example, Δδ¹³C_DIC is always below 0.35 ‰ for Δland, below 450 GtC, and always below 0.23 ‰ for negative Δland estimates. Thus, not only do very few process combinations yield a negative Δland in our 500 000-member ensemble (Fig. 5), but these solutions are also biased low in Δδ¹³C_DIC. Further, solutions with negative Δland also yield ΔCO₂ and $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ values at the very low end of the observational constraints (Fig. 6b, c). Taken together, our framework strongly suggests that the land biosphere sequestered carbon over the deglaciation and that the sequestered amount is likely larger than 450 GtC.

3.3 Carbon and δ¹³C changes in LGM to PI sensitivity simulations

This section is intended for readers interested in a more detailed quantification of deglacial carbon cycle mechanisms. We recall that the model is initialized with a spin-up for PI conditions and simulations started 40 ka under LGM forcing (Fig. 3). Forcing the model from PI into an LGM state by applying standard forcings leads to a shoaling and slight slowdown of the Atlantic meridional overturning circulation (AMOC) of about 4.4 Sv (PI steady state: 18 Sv) and a global cooling of the ocean and the atmosphere of about 1.4 and 3 ^∘C over the glacial, respectively. New production of organic matter decreases from 11.98 GtC yr⁻¹ to a mean value over the glacial of 11.25 GtC yr⁻¹. The oceanic carbon inventory increases at the expense of the atmosphere, yielding a higher ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ concentration during the glacial. Turning to the deglaciation, CO₂ increases by 27.8 ppm from LGM to PI, while δ¹³C_atm and δ¹³C_DIC show little change between LGM and PI (Fig. 7). By applying the standard transient forcings alone, neither the Δδ¹³C constraints in the atmosphere and ocean nor the ΔCO₂ constraint is met, calling for the consideration of additional processes.

Figure 7 Sensitivity of (a) Δδ¹³C_DIC, (b) Δδ¹³C_atm, (c) ΔCO₂, and (d) $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ to changes in Southern Ocean wind stress (SO wind), Southern Ocean gas transfer rate (SO kgas), export rain ratio, shallow water carbonate deposition (coral), remineralization profile (rem. profile), organic weathering flux (org. wea.), and land carbon uptake (land). Sensitivities are shown relative to the standard forcings (white dots). Values in panel (c) correspond to the entries in Table 2. Black crosses with error bars show estimates based on measurements (ΔCO₂: Petit et al., 1999; Monnin et al., 2001; Siegenthaler et al., 2005; Lüthi et al., 2008; Köhler et al., 2017; Δδ¹³C_atm: Schmitt et al., 2012; $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$: Yu et al., 2013; Qin et al., 2017; Luo et al., 2018; Δδ¹³C_DIC: Peterson et al., 2014). Dots in dark grey shading show absolute and in light grey shading relative (to standard forcings) values. Δ indicates the PI–LGM difference.

3.3.1 Sensitivity of carbon cycle processes to changes in additional deglacial processes

The mean δ¹³C signature of the ocean–atmosphere system changes in response to potential deglacial processes (Fig. 7a–b). For simplicity, we focus on the change in mean ocean δ¹³C_DIC (Fig. 7a) as the ocean holds about 20 times more carbon than the atmosphere and interactive land biosphere together. The factorial simulations reveal a Δδ¹³C_DIC (PI–LGM) between −0.3 ‰ and more than 0.6 ‰ when halving or doubling the weathering flux of organic material, respectively. Prescribed changes in Southern Ocean (SO) wind stress inducing changes in ocean circulation, the organic matter remineralization profile, and in the rain ratio lead to a change in the ocean signature of about −0.2 ‰, while changes in the SO air–sea gas transfer rate and in coral reef growth have a small influence on Δδ¹³C_DIC. The influence of these mechanisms on δ¹³C_DIC cannot be ignored when estimating Δland based on reconstructed changes in δ¹³C.

Reconstructions of the change in δ¹³C in the atmosphere and ocean differ in magnitude. The change in the atmosphere (0.1±0.05 ‰) is about one-third of the oceanic change (0.34±0.19 ‰). Changes in the organic weathering flux and terrestrial carbon storage affect Δδ¹³C_DIC and Δδ¹³C_atm in an identical way. Changes in SO wind stress, rain ratio, coral reef growth, and the remineralization of organic material have a similar influence on the two isotope variables. In contrast, changes in the SO air–sea gas transfer rate alter Δδ¹³C_atm, while hardly changing Δδ¹³C_DIC. The strong sensitivity of δ¹³C_atm to changes in SO gas transfer rate is due to the large disequilibrium between the atmosphere and surface ocean. This makes this process, in addition to temperature-driven changes in fractionation, potentially important in setting atmospheric δ¹³C independently from mean ocean δ¹³C_DIC change (see also Fig. C2b).

Considering ΔCO₂ and $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$, the changes in investigated mechanisms significantly influence atmospheric CO₂ and deep Pacific carbonate ion concentration (with the exception of air–sea gas transfer rate). A general relationship becomes apparent. A positive change in ΔCO₂ is always accompanied by a negative change in $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ and vice versa (Fig. 7c–d). Thus, the explanation for the observed deglacial increase in CO₂ of ∼90 ppm requires a deglacial decrease in deep Pacific carbonate ion concentration (Fig. 6d), at least in our model and considering the mechanisms implemented in the standard and factorial simulations. Changes in alkalinity and DIC (PI–LGM) for the standard forcing amount to about 500 Gt and 365 GtC, respectively. In the case of rain ratio, coral reef growth, and changes in the remineralization profile, changes in alkalinity and DIC are largest and occur close to a 2:1 ratio. Changes in land carbon uptake remove carbon from the system, leading to changes in alkalinity via CaCO₃ compensation, whereas in the case of the organic weathering flux, not only carbon is added and/or removed from the system but also nutrients and alkalinity, following Redfield ratios (see Appendix A). Changes in alkalinity and DIC resulting from changes in Southern Ocean wind stress and gas transfer rate are small.

Figure 8 ΔCO₂ for defined sensitivity experiments plotted against (a) Δδ¹³C_atm, (b) Δδ¹³C_DIC, (c) $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$, and (d) Δδ¹³C_atm vs. Δδ¹³C_DIC. Black crosses indicate results with standard forcings alone, and dark grey shading shows ranges for data-based constraints (ΔCO₂: Petit et al., 1999; Monnin et al., 2001; Siegenthaler et al., 2005; Lüthi et al., 2008; Köhler et al., 2017; Δδ¹³C_atm: Schmitt et al., 2012; $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$: Yu et al., 2013; Qin et al., 2017; Luo et al., 2018; Δδ¹³C_DIC: Peterson et al., 2014). For the direction of change in the forcings see Fig. 7. Note that in (a) changes in Δland that yield a positive ΔCO₂ are overlaid by the organic weathering flux line.

3.3.2 Response relationships between different carbon cycle properties

The results for the target variables are plotted against each other to gain further insight into the relative importance of the different mechanisms to meet the observational constraints (Fig. 8). This yields a characteristic slope for each target pair and mechanism. For example, this slope is negative and equates to about −90 ppm ‰⁻¹ for the pair ΔCO₂ and Δδ¹³C_DIC and the mechanism of land carbon storage (grey line in Fig. 8a). Roughly similar slopes are found for variations of the organic weathering flux, SO wind stress, and SO air–sea gas transfer rate. Figure 8a reveals that variations in these four processes cannot prompt the model results to agree with these two observational constraints. The slopes are too small in magnitude and either ΔCO₂ or Δδ¹³C_DIC remains outside the observational range. The four processes are effective in varying Δδ¹³C_DIC but relatively ineffective in varying ΔCO₂. Processes that have a much larger impact on ΔCO₂ than on Δδ¹³C_DIC are related to the CaCO₃ cycle (coral reef growth and rain ratio) and changes in the upper ocean remineralization of particulate organic matter. This implies that variations in Δland, organic weathering, and ocean circulation (SO wind) are required to meet the Δδ¹³C_DIC target, while variations in coral reef regrowth, rain ratio, and upper ocean remineralization depth are necessary to meet the ΔCO₂ target. Indeed, the median changes for the parameters in the observationally constrained ensemble (Fig. C2) imply a deeper remineralization of organic matter and a reduced rain ratio at LGM than PI and a 3 to 4 times higher amount of CaCO₃ deposition during coral reef regrowth than applied in the standard run. These, together with a reduced SO wind stress at the LGM, all contribute to a deglacial increase in CO₂. The associated decrease in Δδ¹³C_DIC from these changes is in most Monte Carlo realizations more than compensated for by a deglacial increase in land carbon inventory and in many realizations by an increased glacial organic weathering flux (Figs. C2, 8f–g). In general, the smaller the change in Δland, the larger the increase required in the glacial flux of organic weathering.

Slopes for the target pair Δδ¹³C_atm−ΔCO₂ resemble the Δδ¹³C_DIC−ΔCO₂ pair (Fig. 8a–b). Variations in SO air–sea gas transfer rate, SO wind stress, Δland, and organic weathering flux show relatively high sensitivity for Δδ¹³C_atm and low sensitivity for ΔCO₂, whereas coral reef regrowth, rain ratio, and upper ocean remineralization changes exert a large impact on ΔCO₂ and affect Δδ¹³C_atm little (Fig. 8b). In contrast to the Δδ¹³C_DIC−ΔCO₂ target pair (Fig. 8a), the slopes for the Δδ¹³C_atm−ΔCO₂ pair do not fall on roughly two characteristic slopes but vary across processes (Fig. 8b).

The clear relationship between oceanic alkalinity and oceanic carbon uptake (Fig. 7c–d) is also visible in Fig. 8c. Higher alkalinity during the glacial (negative $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$) generally results in a larger drawdown of atmospheric CO₂ (positive ΔCO₂). An interesting exception is that a modest deepening of the remineralization of organic matter in the upper ocean leads to large changes in atmospheric CO₂, while at the same time hardly changing deep Pacific ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$, making it a potentially important process to fulfill the ΔCO₂ target without shifting $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ away from the observationally constrained range (Fig. 8c). Accordingly, the remineralization depth of organic matter is deeper at the LGM than PI in all Monte Carlo realizations (Fig. C2e).

Next, we consider the Δδ¹³C${}_{\mathrm{DIC}}-{\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{atm}}$ target pair (Fig. 8d). The difficulty in simultaneously fulfilling the atmospheric and oceanic δ¹³C constraints with a negative Δland is becoming obvious. With negative Δland, Δδ¹³C_atm and Δδ¹³C_DIC are both moved towards negative values relative to standard forcings. Therefore, the isotopic constraints could only be fulfilled with an increased glacial organic weathering flux such that Δδ¹³C_atm and Δδ¹³C_DIC are both moved towards positive values relative to the standard forcings (Fig. 8d). Considering the ΔCO₂ and $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ constraints (Fig. 8c) implies that negative Δland can only be offset by increased glacial organic weathering or rain ratio. However, reaching a ΔCO₂ of 80–100 ppm remains difficult for negative Δland when combining information from all panels of Fig. 8. This explains why only very few solutions with negative Δland and no solutions with Δland smaller than −220 GtC are found in the Monte Carlo simulations (Sect. 3.2).

3.3.3 Sensitivity of the spatial distribution of Δδ¹³C_DIC to changes in deglacial processes

In this section, we describe the spatial distribution of changes in Δδ¹³C_DIC. The reference run with standard LGM forcings (Fig. 9a) shows substantial spatial changes in Δδ¹³C_DIC from internal reorganization linked to ocean circulation changes. Reduced glacial ocean ventilation as evidenced in changes in ideal age leads to the accumulation of light carbon from marine export productivity in the ocean interior. Circulation changes hence explain the variability in the spatial pattern of Δδ¹³C_DIC in Fig. 9a. In particular, the deepening and strengthening of the AMOC during the deglaciation lead to positive Δδ¹³C_DIC changes in the deep Atlantic and to negative Δδ¹³C_DIC changes in the upper Atlantic, consistent with proxy reconstructions.

Figure 9 Simulated deglacial change in the isotopic signature of DIC, Δδ¹³C_DIC (PI minus LGM). Values are displayed along sections across the Atlantic, Southern Ocean, and Pacific and in response to standard deglacial forcings (b–h) due to changes in individual mechanisms as inferred from factorial simulations. Δδ¹³C_DIC is shown for a deglacial carbon uptake of 445 GtC by the land biosphere (b), Southern Ocean wind stress (c), and gas transfer rate (d) scaled by a factor of 0.4 during the glacial relative to PI, rain ratio reduced to 0.045 during the glacial (e), deglacial shallow water carbonate deposition scaled by a factor of 4 (f), a change in the organic matter remineralization between a linear profile during the glacial and the standard depth scaling (g), and for a twofold increase in the organic weathering flux during the glacial compared to PI (h). The absolute change in Δδ¹³C_DIC is shown in panel (a) for the standard forcings (bright colors), while panels (b) to (h) show the difference relative to the standard run (pastel colored).

The sensitivity simulations may be grouped by processes with small and large changes in the spatial pattern of Δδ¹³C_DIC. Δland, coral reef regrowth, and organic weathering flux affect the input of carbon and carbon isotopes into the atmosphere–ocean system, leading to relatively uniform patterns of change (Fig. 9b, f, and h). Changes in SO wind stress, SO gas transfer rate, rain ratio, and upper ocean remineralization profile, on the other hand, mainly lead to a redistribution of carbon and carbon isotopes in the atmosphere–ocean system amplified by sediment interactions (Fig. 9c, d, e, g).

The uptake of carbon by the land biosphere leaves the atmosphere and ocean relatively enriched in δ¹³C (Fig. 9b). As a result of the carbon removal from the coupled system, changes in bottom water concentrations of ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ and O₂ may influence sedimentation fluxes of POC and CaCO₃ locally and lead to small local patterns (see, for instance, negative anomalies in Δδ¹³C_DIC in the North Atlantic deep water (NADW) region of Fig. 9b). Generally, though, a smooth pattern emerges for changes in Δland. Similarly, the doubling of the organic weathering flux during the glacial results in a strong uniform change in Δδ¹³C_DIC (Figs. 9h and 7a). Increased rates of coral reef regrowth have little impact on Δδ¹³C_DIC as the isotopic signature is similar for DIC and deposited CaCO₃ (Fig. 9f). The slight increase in Δδ¹³C_DIC in response to enhanced coral reef regrowth has been attributed to changes in fractionation during marine photosynthesis (Freeman and Hayes, 1992; Menviel and Joos, 2012).

Changes in the other four mechanisms affect the spatial pattern of Δδ¹³C_DIC significantly (Fig. 9c, d, e, g). We note that these patterns result from adjustment processes acting on different timescales. We recall the forcing history for changes in Southern Ocean wind stress, Southern Ocean gas transfer rate, the rain ratio, the remineralization profile, and the organic weathering rate; PI values are set to LGM values at 40 ka, kept constant from 40 to 18 ka, scaled back to the PI value over the termination (18 to 11 ka), and kept constant thereafter (see also Fig. 3 and Table 2). Results are therefore not directly comparable to the step-change experiments described by Tschumi et al. (2011). As discussed in detail by Tschumi et al. (2011) and Menviel et al. (2015), a more poorly ventilated ocean due to reduced wind stress leads to an increase in δ¹³C_atm. Changes in temperature and circulation in response to decreased Southern Ocean wind stress at the beginning of our simulation (Fig. 3f) lead to lower marine oxygen concentrations and an increase in POC sedimentation. Subsequent removal of ¹³C-depleted carbon from the ocean increases δ¹³C_DIC. Restoring wind stress to its initial value over the deglacial reverses the changes. As the sediment feedbacks act on longer timescales, propagation of the signal to the whole ocean is not yet achieved. A complex interplay of processes, such as changes in circulation and subsequent changes in export fluxes, oxygen concentrations, and remineralization of POC, as well as changes in the lysocline and CaCO₃ cycling spatially, overlay each other, yielding the Δδ¹³C_DIC distribution seen in Fig. 9C.

Ocean–sediment interactions also play an important role for δ¹³C_DIC in the case of changes in the Southern Ocean gas transfer rate. Decreasing the gas transfer rate at the beginning of the 40 kyr simulation (Fig. 3f) leads to a strong increase in δ¹³C_atm. The ocean–sediment interactions are very similar to those described for SO wind stress changes above. The sedimentation fluxes of POC and CaCO₃ increase in response to reduced oxygen concentrations from a lower gas transfer rate and changes in the ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ concentration, as does the δ¹³C signature of the flux, leading to a removal of ¹³C from the ocean. From LGM to PI the opposite process takes places (increasing gas transfer back to its initial value); however, due to the long timescales of ocean–sediment interactions, the signal has not propagated to the whole ocean, and in the Atlantic and Pacific negative Δδ¹³C_DIC still prevails (Fig. 9d).

As discussed in detail in Tschumi et al. (2011), responses in δ¹³C from reducing the export rain ratio are small in both the ocean and atmosphere (Fig. 7a–b) and arise from weathering–sedimentation imbalances and changes in the isotopic value of CaCO₃ and POC. This negative perturbation in δ¹³C is not removed entirely over the deglacial and Holocene due to the long ocean–sediment response timescale, leaving the generally negative Δδ¹³C_DIC as seen in Fig. 9e.

Finally, a change to a linear remineralization profile in the upper 2 km during the glacial leads to a dipole pattern in δ¹³C_DIC with enriched values in the upper 1 km of the ocean and depleted values underneath (not shown). Across the deglacial, the process is reversed, yielding the pattern in Δδ¹³C_DIC seen in Fig. 9g from the overlap of different response timescales. Applied changes in the remineralization profile lead to loss of carbon from the atmosphere–ocean system from LGM to PI.

In summary, the carbon and carbon isotope balance in our model simulations can change as a result of altered input fluxes, changes in the surface ocean δ¹³C signatures, which are then reflected in POC and CaCO₃ export fluxes and coral reef regrowth, and altered sedimentation fluxes as a result of altered bottom water ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ and O₂ concentrations. For all seven mechanisms, weathering–sedimentation imbalances and feedbacks exert an impact on the carbon and carbon isotope budget and need to be taken into account on glacial–interglacial timescales. Also, the overlap of different response timescales for atmosphere–ocean and ocean–sediment interactions shapes the patterns in Δδ¹³C_DIC.

4 Discussion

We show that carbon exchange with ocean sediments and the lithosphere biases earlier δ¹³C-based estimates of deglacial change in land carbon storage towards lower values. This finding appears very robust considering modern estimates and reconstructions of ocean–sediment and weathering fluxes and our understanding of the marine carbonate chemistry. Many earlier estimates of land carbon storage change rely on the budget of the stable carbon isotope but neglected the isotopic exchange fluxes with marine sediments and the lithosphere. This approach became popular, perhaps as it is simple and transparent to solve two global budget equations. However, results from idealized pulse–release experiments as well as from transient deglacial simulations show that this neglect is not justified. Several processes that were potentially important for the reorganization of the deglacial carbon cycle influence the mean ocean δ¹³C signature significantly. Their influence cannot be ignored when analyzing the budget of δ¹³C on deglacial timescales. The importance of ocean sediment and lithosphere fluxes in regulating past atmospheric CO₂ and carbon, isotopes, nutrient, and alkalinity inventories in the ocean is also emphasized in previous modeling studies (e.g., Huybers and Langmuir, 2009; Tschumi et al., 2011; Roth and Joos, 2012; Roth et al., 2014; Wallmann et al., 2015; Heinze et al., 2016). Data-based reconstructions of burial fluxes of calcium carbonate and organic matter (Cartapanis et al., 2016, 2018) reveal that these fluxes are more dynamic than previously assumed.

4.1 Uncertainties of this study

There are a number of uncertainties in our approach. We rely on characteristic forcing–response relationship as obtained with the Bern3D model, and these may be different than in reality. The current crop of models, including the Bern3D model, is not able to freely simulate reconstructed variations in atmospheric CO₂ and other biogeochemical variables over the deglacial period. Nevertheless, the Bern3D model simulates the modern distribution of a range of water mass, ventilation, and biogeochemical tracers in good agreement with observational data (Roth et al., 2014). Global inventories and spatial distributions of biogenic CaCO₃, organic carbon, and opal in the sediments, as well as their fluxes between the ocean, marine sediments, and the lithosphere are in agreement with preindustrial data-based estimates (see Appendix B and Table 1). The Bern3D model simulates a weaker and shallower Atlantic meridional overturning circulation at the LGM compared to the PI, and the simulated anomalies in δ¹³C of DIC agree with corresponding reconstructions (see Fig. 9 in Menviel et al., 2012). Regarding dissolved oxygen, qualitative reconstructions show a general decrease in oxygenation at intermediate depths and increase in the deep ocean (Jaccard and Galbraith, 2012; Jaccard et al., 2014; Galbraith and Jaccard, 2015) across the deglaciation. This pattern is also reproduced (not shown) in Bern3D model simulations.

We selected a set of seven generic, archetypical processes that are varied in deglacial simulations in addition to the processes explicitly implemented in the Bern3D model. These archetypical processes are intended to represent the space of multi-proxy response relationships (Fig. 8) for all the different processes that plausibly influenced deglacial carbon cycle changes. This is a simplification, though we argue that the seven selected processes roughly cover the plausible range of multi-proxy responses (Fig. 8). In addition, uncertainties in the importance of these processes are explicitly considered in our Monte Carlo approach and contribute to the uncertainty range in our estimate of Δland.

For example, there are various carbon reservoirs with an isotopic light signature very similar to that of plant and soil carbon. These include organic carbon stored in ocean sediments and shelves (Broecker, 1982; Ushie and Matsumoto, 2012; Cartapanis et al., 2016), dissolved organic carbon (DOC) in the ocean (Hansell, 2013), and gaseous CO₂ stored in the unsaturated zones of aquifers from the decomposition of organic material (Baldini et al., 2018). Any carbon release or uptake from these reservoirs might be wrongly attributed to a change in the carbon inventory of the land biosphere given their very similar δ¹³C signature and similar multi-proxy response relationships. Indeed, changes in these reservoirs have been invoked as an alternative explanation for the reconstructed deglacial change in marine δ¹³C (Broecker, 1982; Zimov et al., 2006, 2009; Zech et al., 2011; Kemppinen et al., 2018). The deglacial change in CO₂ stored in the unsaturated zone of aquifers is in the range of 12 to 86 GtC (Baldini et al., 2018) and is small compared to the estimated change in land carbon of around 850 GtC. Perhaps more important, the modern inventory of refractory DOC is about 650 GtC (Hansell, 2013) and assumed constant in Bern3D, while changes in the small inventory of labile DOC are explicitly simulated. It is unclear how the DOC pool varied over the past. Further, Cartapanis et al. (2016) suggest an approximately linear decrease in the flux of organic carbon transferred to deep ocean sediments from around 28 GtC kyr⁻¹ at the LGM to around 17 GtC kyr⁻¹ during the Holocene, corresponding to a cumulative imbalance of around 55 GtC. The change in organic carbon stocks stored in sediments on shelves is uncertain and may amount to several hundred gigatons of carbon over the deglaciation. Here, we varied organic-like carbon input by weathering in the unconstrained ensemble over a range corresponding to a cumulative PI–LGM anomaly of +780 to −1560 GtC. These anomalies in input, though not directly comparable to inventory changes, are large in comparison to the inventory changes discussed here. In summary, a potential misattribution of changes in other organic carbon reservoirs to changes in Δland is included in our uncertainty estimate of Δland.

Our approach is limited in that we only consider the change between LGM and PI for four relevant proxy variables. This restriction allowed us to build a cost-efficient and nonlinear emulator to explore a very large range of parameter combinations. Future efforts to refine estimates of Δland may explicitly consider the spatial and temporal evolution for a large number of proxies.

4.2 Contribution of deglacial mechanisms to the CO₂ rise

Turning to the role of different mechanisms for the deglacial increase in atmospheric CO₂, we find that ocean circulation changes only partly explain the deglacial CO₂ rise. Ocean circulation changes are not able to simultaneously explain glacial–interglacial changes in δ¹³C_DIC and CO₂ in our model (see Fig. 8a), in agreement with earlier studies (e.g., (Tschumi et al., 2011; Heinze et al., 2016)).

A modest and plausible (Bendtsen et al., 2002; Matsumoto, 2007) deepening of the remineralization of organic matter in the upper ocean (see Fig. C2e) leads to large changes in atmospheric CO₂, while at the same time hardly changing deep Pacific ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$. This also holds for other biological mechanisms such as an increased nutrient utilization in response to iron fertilization or an elevated phosphorus input to the ocean (Menviel et al., 2012). This renders this class of processes potentially important to explain part of the deglacial CO₂ increase, as reconstructions suggest small LGM to PI changes in deep Pacific ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ concentrations (Yu et al., 2014). Accordingly, in all Monte Carlo realizations the remineralization depth of organic matter is deeper at the LGM than PI in our approach.

Our solutions also point to a significant role of alkalinity-based mechanisms for the deglacial CO₂ increase (see Fig. C2d). The constrained ensemble yields a large extra burial of CaCO₃ as well as an increase in the rain ratio of CaCO₃ to organic matter by particle export over the deglaciation. This is qualitatively consistent with the reconstruction of Cartapanis et al. (2018), who suggest that the CaCO₃ burial flux (below 200 m) increased by about 0.2 GtC yr⁻¹ from the LGM to the Holocene.

4.3 Comparison to other estimates of Δland

Δland estimates range from −400 GtC, e.g., a larger terrestrial carbon inventory during the glacial (Zimov et al., 2006, 2009; Zech et al., 2011; Kemppinen et al., 2018), to 1500 GtC (Adams and Faure, 1998) based on a variety of methods. The low end is based on reconstructions of the carbon content in paleosols (Zimov et al., 2006, 2009; Zech et al., 2011). These reconstructions are scarce, of local nature, and restricted to cold high-altitude or high-latitude regions, rendering extrapolations to the global scale speculative. In a recent data synthesis, Lindgren et al. (2018) suggest that the loss of carbon from Northern Hemisphere permafrost from the LGM to PI was lower than initially suggested by Zimov et al. (2009) and that total land carbon storage in northern middle and high latitudes increased by about 400 GtC over the deglaciation. Kemppinen et al. (2018) suggest a negative Δland based on results from an Earth system model of intermediate complexity. However, these authors did not consider any isotopic constraints. Our results suggest negative Δland values to be highly unlikely and Δland to be likely larger than 450 GtC. Otherwise, solutions for δ¹³C_DIC are significantly biased low compared to the data-based reconstruction.

The high end of the Δland range, implying much smaller land carbon storage during the LGM than today, is based on estimates relying on pollen records (Adams et al., 1990; Van Campo et al., 1993; Crowley, 1995; Adams and Faure, 1998). Pollen and macrofossils recorded in various archives hold information on past vegetation distribution but do not allow for constraining carbon inventories in soils or in vegetation, rendering the derivation of Δland estimates based on pollen data difficult (Hoogakker et al., 2016). Ice core CO₂ and δ¹³C_atm records constrain the preindustrial terrestrial carbon uptake over the Holocene to about 250 GtC (Elsig et al., 2009). Assuming the total estimate of 1500 GtC by Adams and Faure (1998) was right, this would leave more than 1200 GtC to be transferred to the terrestrial biosphere between the LGM and the beginning of the Holocene. Considering possible processes that could achieve such a large transfer seems difficult. Generally this might pose an upper limit to Δland.

Our estimate of Δland appears consistent with a recent biome-based reconstruction of soil carbon storage (Lindgren et al., 2018). Lindgren et al. (2018) report an increase in soil carbon inventory of about 400 GtC for the region north of about 30^∘ N and considering changes in permafrost, mineral, and peatland, as well as in loess and subglacial soil. This is in agreement with our Δland estimate when assuming an increase in carbon stored in vegetation and in tropical and Southern Hemisphere ecosystems and soils. Furthermore, modeling studies have generally yielded positive Δland estimates of about 400–900 GtC (Joos et al., 2004; Prentice et al., 2011; Brovkin et al., 2012; O'Ishi and Abe-Ouchi, 2013; Ganopolski and Brovkin, 2017; Davies-Barnard et al., 2017).

In summary, many lines of evidence point to a positive Δland, and we suggest a likely range of about 450 to 1250 GtC. Yet, the uncertainty in this estimate remains large owing to both uncertainties in processes and in observational constraints. Unfortunately, the necessity of including exchange with sediments and the lithosphere and other deglacial processes in the isotopic budget adds complexity. This tends to increase the uncertainty range compared to simplified and biased assessments that consider the isotopic budget in the closed ocean–land–atmosphere system.

5 Summary and conclusions

We used a Bayesian approach to constrain the change in the land biosphere carbon inventory between the Last Glacial Maximum and the preindustrial period. Four data-based constraints are applied: the deglacial change in the δ¹³C signature of atmospheric CO₂ and of dissolved inorganic carbon in the ocean; and the deglacial change in atmospheric CO₂ and in deep Pacific carbonate ion concentrations. The strength of generic, archetypical mechanisms for deglacial biogeochemical changes was varied in 500 000 simulations in an emulator built from a large suite of Bern3D model simulations. Carbon, nutrient, alkalinity, and carbon isotope exchange with the lithosphere and ocean sediments is explicitly taken into account, in contrast to earlier studies that applied the isotopic budget of ¹³C to estimate deglacial change in land biosphere carbon.

We demonstrate in idealized and transient deglacial simulations that isotopic exchange with ocean sediments and the lithosphere is important for the budget of ¹³C on timescales of the deglaciation. We find that the carbon stocks in the land biosphere were likely around 450 to 1250 GtC (±1 standard deviation) larger at preindustrial times than at the Last Glacial Maximum. The median estimate is ∼850 GtC. This is much larger than earlier estimates based on the budget of the stable carbon isotope ¹³C. These earlier studies neglected interactions with the lithosphere and ocean sediments and the influence of other deglacial carbon cycle processes on the isotopic budget. This neglect biases their estimate significantly low. Our multi-proxy approach suggests that a combination of different mechanisms contributed to the deglacial increase in CO₂. These include, aside from ocean circulation, temperature and salinity changes, a net removal of alkalinity from the surface ocean through increased burial of calcium carbonate in coral reefs and ocean sediments, and potentially increasing export of calcium carbonate from the surface to the deep ocean. Further, changes in biological mechanisms such as a deglacial shoaling of the remineralization depth for organic matter may have been instrumental in increasing atmospheric CO₂. The results demonstrate that ocean sediments and the weathering–burial cycle are an integral part of the Earth system playing a fundamental role on glacial–interglacial timescales.

Appendix A: The Bern3D model and pulse experiments

A1 The Bern3D model

The Bern3D EMIC features a three-dimensional geostrophic ocean (Müller et al., 2006; Edwards et al., 1998) with an isopycnal diffusion scheme and Gent–McWilliams parameterization for eddy-induced transport (Griffies, 1998). The model includes a thermodynamic sea-ice component coupled to a single-layer energy–moisture balance atmosphere (Ritz et al., 2011). Further, a sediment module (Tschumi et al., 2011; Heinze et al., 1999) and a four-box terrestrial biosphere (Siegenthaler and Oeschger, 1987) are coupled to the model.

The horizontal resolution is 41×40 grid cells (Roth et al., 2014; Battaglia and Joos, 2018a, b) and is the same for the ocean, atmosphere, and sea-ice components. In the vertical, the ocean has 32 logarithmically spaced layers. Wind stress at the surface is prescribed following the NCEP/NCAR monthly wind stress climatology (Kalnay et al., 1996). Carbonate chemistry and the air–sea gas exchange for CO₂ and ¹⁴CO₂ are implemented according to OCMIP-2 protocols (Najjar and Orr, 1999; Orr et al., 1999) with updates for the ¹⁴C standard ratio and half-life (Orr et al., 2017), the calculation of the Schmidt number (Wanninkhof, 2014), and the carbonate chemistry (Orr and Epitalon, 2015). In order to match observational estimates of natural and bomb-produced ¹⁴C, the global mean air–sea transfer is reduced by 19 % compared to OCMIP-2 in order to match observation-based radiocarbon estimates (Müller et al., 2008). Another update with respect to the OCMIP-2 protocols concerns the gas transfer velocity, which scales linearly with wind speed following Krakauer et al. (2006). Analog and consistent formulations are used for air–sea exchange of oxygen and ¹³CO₂.

Biogeochemical cycling in the model is detailed in Tschumi et al. (2011) and Parekh et al. (2008) with further documentation of results in follow-up studies (e.g., Menviel et al., 2012, 2015; Menviel and Joos, 2012; Roth and Joos, 2012; Roth et al., 2014; Battaglia et al., 2016; Battaglia and Joos, 2018b, a). Dissolved inorganic carbon and semi-labile organic carbon (DIC, DOC), the corresponding isotopic forms (DI¹³C, DO¹³C, DI¹⁴C, DO¹⁴C), and alkalinity (Alk), phosphate (PO₄), oxygen (O₂), iron (Fe), silica (Si), and an ideal age tracer are explicitly transported by advection, diffusion, and convection. New production of organic matter is limited to the euphotic zone in the uppermost 75 m and calculated as a function of light availability, temperature, and phosphate and iron availability (Doney et al., 2006). Bacteria and plankton are not explicitly represented. Two-thirds of the new production is transferred to the pool of dissolved organic carbon (DOC) and the remainder is exported as particulate organic carbon (POC) (Tschumi et al., 2008). The biological fluxes of C, PO₄, Alk, and O₂ and the corresponding elemental ratios in DOC and POC are linked by fixed Redfield ratios ($\mathrm{P}:\mathrm{Alk}:\mathrm{C}:{\mathrm{O}}_{\mathrm{2}}:\mathrm{Alk}=\mathrm{1}:\mathrm{17}:\mathrm{117}:-\mathrm{170}$; Sarmiento and Gruber, 2006). DOC decays with a lifetime of 0.5 yr. POC is remineralized following a Martin's curve (Martin et al., 1987) given by

$$\begin{array}{}\text{(A1)}& F_{\mathrm{POC}}\left(z\right)=F_{\mathrm{POC}}\left(z_{\mathrm{0}}\right)\times {\left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)}^{-\mathit{\alpha }}\text{for}z>z_{\mathrm{0}},\end{array}$$

where z₀ refers to the reference depth of 75 m (see also Roth et al., 2014). Export of calcium carbonate (CaCO₃) and of opal is computed from new production and availability of dissolved silica in the euphotic zone. The export rain ratio of CaCO₃ to POC is set to 0.083 in the absence of dissolved silica, while the production ratio of CaCO₃ to POC is reduced at the favor of opal production in regions with abundant availability of silicic acid. Particulate CaCO₃ and opal are remineralized below the euphotic zone with an e-folding depth scale of 5066 and 10 000 m, respectively.

A 10-layer sediment diagenesis module (Tschumi et al., 2011; Heinze et al., 1999) is used to explicitly calculate transfer fluxes to and redissolution fluxes from reactive sediment to the water column as well as loss fluxes to the lithosphere for nutrients, carbon, and carbon isotopes. In equilibrium, loss fluxes to the lithosphere are equal to input fluxes by weathering. The sediment module dynamically calculates the transport, redissolution, remineralization, and bioturbation of solid material, the pore-water chemistry, and diffusion in the top 10 cm of the sediment. Four solid tracers (CaCO₃, opal, POC, clay) and seven tracers in the pore water (DIC, DI¹³C, DI¹⁴C, alkalinity, phosphate, oxygen, and silicic acid) are modeled. The dissolution of CaCO₃ depends on the pore-water ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ concentration. The oxidation rate of POC, on the other hand, depends on oxygen concentrations in the pore water and the weight fraction of POC in the solid phase of the sediment. Denitrification is not considered in this model version, but O₂ is not consumed below a threshold, somewhat reflecting the process of denitrification without modeling ${\mathrm{NO}}_{\mathrm{3}}^{-}$. The respective reaction rate parameters for CaCO₃ dissolution and POC oxidation are global constants (see Roth et al., 2014). The model assumes conservation of volume; i.e., the entire column of the sediments is pushed downwards if deposition exceeds redissolution into pore waters. Any solid material that is pushed out of the diagenetic zone (top 10 cm) disappears into the subjacent diagenetically consolidated zone (burial or loss flux) (see Tschumi et al., 2011, for more details). Carbonate chemistry within sediment pore waters is calculated as in the ocean by using the MOCSY routine of Orr and Epitalon (2015).

Weathering fluxes of PO₄, Alk, DIC, DI¹³C, and silicic acid are added uniformly to the coastal surface ocean. The global weathering fluxes are set equal to the burial fluxes diagnosed at the end of the model spin-up. In the standard model setup these diagnosed input fluxes are kept constant. The pulse experiments include a setting in which the input fluxes vary as a function of global mean surface air temperature and CO₂ following Colbourn et al. (2013). One of the factorial sensitivity runs assesses the impact of altered input fluxes of organic material across the deglacial. This weathering–burial cycle and the associated burial–nutrient feedbacks (Tschumi et al., 2011) are important for the mass balance of ¹³C and for the removal of marine perturbations in nutrients, carbon, and isotopes on deglacial timescales. In other words, these processes directly affect the estimate of the change in land carbon from the ¹³C mass balance.

The atmosphere exchanges carbon and carbon isotopes with a four-box land biosphere model (Siegenthaler and Oeschger, 1987). It includes a mass of 2220 GtC and represents the carbon that is actively exchanged with the atmosphere. This four-box reservoir is only used to simulate the dilution of an atmospheric isotopic perturbation by the land biosphere but not to address changes in land carbon stocks. The carbon inventory of the land biosphere, and thus the dilution of an isotopic perturbation, is about 20 times smaller than that of the ocean–atmosphere–sediment system.

¹³C stocks and flows are modeled in the atmosphere–ocean–sediment–land biosphere system. ¹³C fluxes to the lithosphere associated with the burial of POC and CaCO₃ and weathering input are explicitly simulated. ¹³C fractionation is considered for air–sea gas transfer, carbonate chemistry, the formation of CaCO₃, POC, DOC, and photosynthesis on land, while no fractionation is assumed during remineralization on land and in the ocean. Gross air-to-sea and gross sea-to-air fluxes of ¹³CO₂ are implemented considering kinetic fractionation (Siegenthaler and Oeschger, 1987) and strongly temperature–dependent equilibrium fractionation between the various carbonate species (Mook, 1986). Isotopic fractionation during CaCO₃ formation is small and computed following Mook (1986). This results in an isotopically heavy signature of around 2.9 ‰ for CaCO₃. Fractionation during photosynthesis in the ocean and thus between dissolved CO₂ ([CO₂]) and POC and DOC is calculated according to Freeman and Hayes (1992). The fractionation increases logarithmically with [CO₂] in surface water and results in an isotopically light signature for POC and DOC of around −20 ‰. A lowering of [CO₂], or correspondingly of pCO₂, by about 10 % yields a change in the isotopic signature of POC by +0.55 ‰. The burial of POC and CaCO₃ results in an isotopic signature of the burial flux of around −9 ‰, intermediate between the isotopically light POC and the isotopically heavy CaCO₃ signature. In equilibrium, the burial flux is compensated for by a weathering input flux of equal amount and signature. On land, a constant fractionation of 18.1 ‰ is applied for simplicity.

A2 Initialization of Bern3D simulations

The model is spun up over 60 000 years to a preindustrial equilibrium corresponding to 1765 CE boundary conditions similar to the approach outlined in Roth et al. (2014). CO₂ is set to 278 ppm and δ¹³C of CO₂ to −6.305 ‰. The loss of tracers due to sedimentary burial is compensated for during the spin-up by variable weathering fluxes in order to conserve oceanic inventories of tracers. After the system has equilibrated, the weathering fluxes are kept constant at the values diagnosed at the end of the spin-up in the standard setup.

A3 Bern3D pulse experiments

Idealized experiments with a pulse-like carbon removal from the atmosphere are conducted in (i) a closed system, (ii) an open system, and (iii) an open system model configuration allowing for weathering feedbacks. In (i), only the atmosphere–ocean–land biosphere model components are used and all tracers are conserved within these three reservoirs. In (ii), the sediment module is added to the configuration outlined above and the carbon and tracer inventories within the atmosphere–ocean–land biosphere system are allowed to vary due to weathering–burial imbalances. In (iii), weathering fluxes vary in response to changes in atmospheric temperature and CO₂ concentrations following the parameterizations given in Colbourn et al. (2013), accounting for changes in CaCO₃ and CaSiO₃ weathering on land.

Simulations are started from the end of the PI spin-up, and atmospheric CO₂ and δ¹³C are computed in a prognostic way based on air–sea and air–land fluxes. In year 100 of the simulation, 100 GtC of carbon with a δ¹³C signature of −24 ‰ is instantaneously removed from the atmosphere, mimicking an immediate regrowth of the terrestrial biosphere. After the pulse, the model is run for 100 kyr and anomalies are expressed relative to a control run. Control simulations show negligible changes in CO₂ and δ¹³C as model drift is small.

The implied change in the land biosphere carbon stock based on a closed system assumption is calculated for all three model settings by solving the following equations for ΔM_TER, i.e., Δland:

$$\begin{array}{}\text{(A2)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }(M_{\mathrm{A}}+M_{\mathrm{B}}+M_{\mathrm{O}})=\mathrm{0},\text{(A3)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }({\mathit{\delta }}_{\mathrm{A}}{M}_{\mathrm{A}}+{\mathit{\delta }}_{\mathrm{B}}{M}_{\mathrm{B}}+{\mathit{\delta }}_{\mathrm{O}}{M}_{\mathrm{O}})=\mathrm{0},\end{array}$$

with M representing the inventory, δ the isotopic signature, Δ the respective time interval, and the subscripts A, B, and O referring to the atmosphere, terrestrial biosphere, and ocean, respectively. The calculated Δland is then compared to the prescribed carbon uptake. This allows us to investigate the validity of the closed system assumption on different timescales.

Appendix B: Sediment tuning

The aim of the sediment tuning was to improve the representation of export, deposition, and dissolution fluxes of POC, CaCO₃, and opal, as well as the distribution of weight fractions of these tracers in the sediment module of the Bern3D model after the change of the horizontal model resolution introduced in Roth et al. (2014). One concern was the total amount of CaCO₃ in the reactive uppermost 10 cm of the sediment represented in the sediment module. Prior to tuning, the global CaCO₃ content in the model was ∼400 GtC. Using observational constraints on modern global bulk mass accumulation rates (MARs) (Cartapanis et al., 2016), global carbonate MAR (Dunne et al., 2012; Jahnke, 1996), and the density, bulk sediment, and porosity in the uppermost 10 cm as used in the model, we estimate a target CaCO₃ stock of ∼1000 GtC. The second concern was the opal burial that used to be ∼3 Tmol Si yr⁻¹ but with new observational data suggesting a burial flux of ∼6.5 Tmol Si yr⁻¹ (Tréguer and De La Rocha, 2013). Further, all other sediment-related fluxes of POC, CaCO₃, and opal were tuned with regard to observational constraints (see Table 1).

Table B1 Old and new parameter values in the Bern3D model that were varied in the tuning of the sediment component.

The tuning was carried out in two steps. In the first step, three parameters were varied: the ratio of Ca:P in calcifiers, the redissolution length scale for the calcite profile, and the reactivity for calcite in the sediment. In the second step, the reactivity of calcite and the diffusion coefficient in the sediment pore water were varied. For both steps, a 100-member LHS ensemble was run and skill scores calculated as outlined in Steinacher et al. (2013). For the calculation of the skill scores, ocean data from GLODAP v.2 (Lauvset et al., 2016) and the World Ocean Atlas 2013 version 2 (Locarnini et al., 2013; Zweng et al., 2013; Garcia et al., 2013a, b) datasets were used. For export, deposition, and burial of POC, CaCO₃, and opal values from Battaglia et al. (2016), Tréguer and De La Rocha (2013), Sarmiento and Gruber (2006), Milliman and Droxler (1996), and Feely et al. (2004) were used. For the sediment data compilations by Cartapanis et al. (2016, 2018) were used to compare the model performance. The distribution of CaCO₃ in the sediments after the tuning compared to observational data is shown in Fig. B1. Export, deposition, and burial fluxes for POC, CaCO₃, and opal and the global sediment inventories in the model with respective observational ranges from the literature are given in Table 1 and the old and new parameter values in Table B1. The distribution of other tracers is comparable to before and as outlined in Roth et al. (2014).

Figure B1 (a) Distribution of CaCO₃ in the Bern3D model after tuning of the sediment component under preindustrial boundary conditions compared to (b) sediment core-top CaCO₃ data (Cartapanis et al., 2018). Core-top data were gridded onto the Bern3D grid and are a mean of the top 10 cm over the past 6 kyr. Grey indicates land (within coastlines) or missing data.

Appendix C: Additional figures

Figure C1 Regression of emulator results for the four data-based constraints against results from a Bern3D LHS ensemble with the same parameter combinations (as outlined in Sect. 2.2). The regression given in the figure is used to correct emulator results for nonlinearities (Eq. 2).

Figure C2 Histograms (colored thin line, normalized) and kernel density estimates (colored thick line) of applied changes in the seven mechanisms (a–g) for all emulator results that fulfill the four observational constraints. The colors of the lines correspond to subsamples of these results according maximum changes in Δland of 1500, 1000, 750, 450, and 0 GtC. The thin red line gives the histogram of the prior distribution and the vertical thick red line the standard parameter value. For details on the mechanisms see Sect. 2.2.