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 ( (cid:49) 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 ﬂux imbalances, and the inﬂuence of the transient deglacial reorganization on the isotopic budgets. We show this simpliﬁcation to signiﬁcantly reduce (cid:49) land in simulations using the Bern3D Earth System Model of Intermediate Complexity v.2.0s. We constrain (cid:49) land to ∼ 850 GtC (median estimate; 450 to 1250 GtC ± 1SD) by using reconstructed changes in atmospheric δ 13 C, marine δ 13 C, deep Paciﬁc carbonate ion concentration, and atmospheric CO 2 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 ﬂuxes involving marine organic matter to explain deglacial change and suggests a major upward revision of earlier isotope-based estimates of (cid:49) land.


Introduction
Atmospheric CO 2 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 2 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 2 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;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 2 variability emerge, these models are not yet able to reproduce variations in important proxy data or to represent the timing of CO 2 changes over the last glacial termination adequately (Brovkin 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.
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 13 C. The argument is that the uptake of isotopically light carbon by the land biosphere causes a corresponding measurable perturbation in the average δ 13 C value in the oceanatmosphere system. The 13 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 δ 13 C sediment records (Duplessy et al., 1984;Curry et al., 1988;Bird et al., 1994Bird et al., , 1996, but recently, more comprehensive δ 13 C data compilations (Oliver et al., 2010;Peterson et al., 2014) have become available. These were also used to determine mean ocean δ 13 C and land by applying a dynamic ocean model (Menviel et al., 2017).
All previous 13 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 2 rise and that large changes in the calcium carbonate (CaCO 3 ) 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 δ 13 C in the coupled atmosphere-surface ocean system to ocean sediments (Tschumi et al., 2011;Roth et al., 2014;Cartapanis et al., 2016). A second reason relates to the chemical buffering of land-induced CO 2 perturbations to marine CaCO 3 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 13 C between the ocean and reactive marine sediments. Third, and of even greater importance with regard to δ 13 C, are changes in the cycling and weathering-sedimentation imbalances of organic carbon that have the potential to exert a large impact on the 13 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 δ 13 C to infer changes in land carbon stocks.
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 δ 13 C DIC and for other proxy targets (PI-LGM differences) to the changes in individual key deglacial carbon cycle mecha-nisms 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 13 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 2 and δ 13 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 wholeocean δ 13 C DIC signature. All previous δ 13 C-based estimates of land neglected the influence of these processes and appear systematically biased towards lower estimates.

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 spinup 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., , 2018Tré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 3 , 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). 13 C stocks and exchanges are modeled in the atmosphereocean-sediment-land biosphere system. 13 C fluxes to the lithosphere associated with the burial of POC and CaCO 3 and weathering input to the ocean are explicitly simulated. The burial of POC and CaCO 3 results in an isotopic signature of the burial flux of around −9 ‰, intermediate between the isotopically light POC (∼ −20 ‰) and the isotopically heavier CaCO 3 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 3 export fluxes from the surface ocean contribute to the burial flux. Figure 2 shows modeled δ 13 C in a section through the Atlantic, Southern Ocean, and Pacific under preindustrial (1765 CE) boundary conditions with a prescribed atmospheric δ 13 C signature of −6.305 ‰. Compared to measured data corrected for the Suess effect (Eide et al., 2017) it is visible that modeled δ 13 C values are biased towards high values by a seemingly constant offset of about 0.4 ‰, while the δ 13 C patterns are faithfully captured by the model.  (Joos and Spahni, 2008) and variations in orbital parameters (Berger, 1978). Radiative forcing of CO 2 is prescribed and biogeochemistry is determined with interactive CO 2 . Ice sheet extent and related changes in albedo are prescribed based on benthic δ 18 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 , 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 13 C from the uppermost ocean grid cells. These forcings  are termed standard forcings.    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 2 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.

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 sen-sitivities 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 Table 1. Export, burial, and deposition fluxes, as well as sediment inventories of CaCO 3 , 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.  Sarmiento and Gruber (2006), d Milliman and Droxler (1996), e Feely et al. (2004).
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 2 versus changes in whole-ocean δ 13 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 2 changes, and a reduced Antarctic overturning circulation has been invoked to, at least partly, explain low glacial atmospheric CO 2 (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;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 (> 48 • 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. Conse-quently, 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 δ 13 CO 2 (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 2 is discussed elsewhere (Sigman and Boyle, 2000;Brovkin et al., 2007;Sigman et al., 2010; and is part of the standard forcing. Oceanic carbonate. The partial pressure of CO 2 in the ocean depends on alkalinity and processes that alter surface ocean alkalinity, thus influencing atmospheric CO 2 . The sedimentary burial of CaCO 3 removes twice as much alkalinity as carbon from the ocean, leading to oceanic CO 2 outgassing. Increased CaCO 3 burial over the deglacial has been put forward as the coral reef hypothesis (Berger, 1982) to explain increases in atmospheric CO 2 . 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 3 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 2 . Here, we scale the reconstructed deposition history based on Vecsei  and Berger (2004) (Fig. 3d) by applying a constant scaling to vary CaCO 3 deposition within the published range.
Potential changes in the rain ratio (CaCO 3 : C org ) affect both alkalinity and CO 2 (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 3 relative to C org from surface waters has similar consequences for CO 2 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 cy- 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.

Mechanisms
Parameter values Southern Ocean wind stress 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 Southern Ocean gas transfer rate 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1. cling of organic matter within the ocean, the whole-ocean nutrient inventories, and thereby surface ocean nutrient and DIC concentrations as well as atmospheric CO 2 . Several studies have suggested changes in the remineralization length scale of organic carbon due to colder ocean temperatures (see, e.g., Matsumoto, 2007;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 Roth et al., 2014), resulting in a decrease in atmospheric CO 2 (Matsumoto, 2007;. Weathering-burial feedbacks lead to an amplification of changes in atmospheric and oceanic δ 13 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  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., , 2018. A negative balance between the riverine input of PO 4 , DIC, DI 13 C, and Alk, originally buried as organic material, and sedimentary burial of particulate organic matter leads to a transient decrease in atmospheric CO 2 and also affects atmospheric and marine δ 13 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 2 . To account for this possibility, we varied the input of PO 4 , DIC, DI 13 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;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 δ 13 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 termina-856 A. Jeltsch-Thömmes et al.: Low terrestrial carbon storage at the Last Glacial Maximum tion (18-11 ka). The δ 13 C signature of terrestrial carbon is set to −24 ‰.

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 2 ( CO 2 ), (ii) atmospheric δ 13 C(CO 2 ) ( δ 13 C atm ), (iii) mean ocean δ 13 C(DIC) ( δ 13 C DIC ), and (iv) deep (equatorial) Pacific carbonate ion concentration ( CO 2− 3 ) ( 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 2 constrains the integrated effect of deglacial carbon cycle changes. We allow for a relatively wide uncertainty range of 20 ppm for CO 2 , which is admittedly wider than the proxy uncertainty, to avoid an overfitting of model outputs. The δ 13 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 δ 13 C DIC of 0.34 ± 0.19 ‰, and we use their range (0.15 ‰ to 0.53 ‰) as a constraint. Proxy reconstructions of CO 2− 3 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 CO 2− 3 in the western tropical Pacific across the last glacial termination. In a recent study, Luo et al. (2018) report, based on CaCO 3 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 CO 2− 3 , 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 CO 2− 3 is chosen here from −6 to 5 µmol kg −1 .

Emulator and Monte Carlo setup
Sensitivities of the four target proxies ( CO 2 , δ 13 C atm , δ 13 C DIC , and CO 2− 3 ) 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 T i ( p i ). These sensitivity functions are used to build a simple representation to compute changes in target variables for 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 δ 13 C and CO 2− 3 , respectively. The assumed uncertainty range in CO 2 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.

Constraint
Uncertainty range References 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: 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 2 , 0.06 ‰ for δ 13 C atm , 0.02 ‰ for δ 13 C DIC , and 4 µmol kg −1 for CO 2− 3 . 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.

Pulse experiments
We start by investigating how atmospheric CO 2 and δ 13 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 (δ 13 C = −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 13 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 2 , followed by a typical impulse-response recovery (e.g., Maier-Reimer and Hasselmann, 1987;Siegenthaler and Joos, 1992). Atmospheric CO 2 recovers as the ocean starts to outgass CO 2 (Fig. 4). In a closed system, atmospheric CO 2 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 13 C. Initially, the δ 13 C atm perturbation is removed much faster than the atmospheric CO 2 perturbation as exchange with the ocean and land biosphere dilutes the imposed isotopic perturbation ( Fig. 4a; Joos and Bruno, 1996). The oceanic δ 13 C increases and equilibrates after about 2 kyr. The difference in the response of CO 2 and δ 13 C is related to their different properties. CO 2 represents a concentration, whereas δ 13 C stands for the isotopic ratio of 13 C to 12 C. The capacity of the ocean to remove an atmospheric perturbation in CO 2 , and equally in 12 CO 2 and 13 CO 2 , 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 13 CO 2 / 12 CO 2 is hardly affected by this chemical buffering; the buffering affects 13 CO 2 in the nominator and 12 CO 2 in the denominator about equally, leading to a near cancellation of its effect on the ratio. The resulting change for the new equilibrium in δ 13 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 2 (see Archer et al., 1998) is similar to the closed system except that CaCO 3 compensation leads to a further recovery of atmospheric CO 2 until equilibrium is approached with an e-folding timescale of about 14 000 years. CO 2 equilibrates at around 3.5 ppm below the initial steady state. For the first thousand years, the evolution of δ 13 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 δ 13 C perturbation in the oceanatmosphere-land biosphere system is increasingly and eventually completely removed by open system processes. Figure 4b shows the δ 13 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 δ 13 C than expected from closed system modeling. The resulting difference in the δ 13 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 13 C budget in the open system changes due to imbalances between input and burial fluxes of POC and CaCO 3 and due to changes in the δ 13 C signatures of the respective fluxes. The uptake of isotopically light land carbon from the atmosphere leads to an increase in the δ 13 C signature of surface waters and through that in the δ 13 C signatures of exported and eventually buried POC and CaCO 3 . The mean δ 13 C signature of the burial flux amounts to −8.7 ‰ (δ 13 C CaCO 3 : 2.9 ‰) and C org (δ 13 C POC : −20.2 ‰) for the first 10 kyr compared to −9.1 ‰ for the weathering input. This difference tends to mitigate the δ 13 C DIC perturbation. Further, carbon-climate feedbacks cause changes in temperature and ocean circulation, which affect the export of POC and CaCO 3 . Burial fluxes of POC further depend on O 2 concentrations and burial fluxes of CaCO 3 on CO 2− 3 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 δ 13 C signature, while changes in the CaCO 3 cycle play a smaller but non-negligible role. Overall, 13 C is lost from the atmosphere-ocean system, and δ 13 C DIC is decreasing.
In the third pulse experiment, changes in global mean air temperature and atmospheric CO 2 concentration affect weathering fluxes of CaCO 3 and CaSiO 3 following Colbourn et al. (2013). The implemented feedbacks cause a reduction in modeled weathering fluxes in response to decreased atmospheric CO 2 and temperature over the course of the experiment. Decreased weathering fluxes of CaCO 3 and CaSiO 3 lead to a decline in alkalinity supplied to the ocean, result- ing in a net transfer of carbon to the atmosphere. This further reduces the initial CO 2 perturbation (Fig. 4a). The mean δ 13 C value related to weathering is −9.1 ‰, thereby leading to a slightly higher δ 13 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 δ 13 C budgets on deglacial timescales. Taken at face value, the results suggest that the classical δ 13 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 glacialinterglacial transition. This reorganization affected open system processes and the 13 C budget in addition to carbon uptake by the land biosphere. In the next section, we estimate land by taking into account deglacial processes.

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 000member 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 δ 13 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 δ 13 C DIC or δ 13 C atm together with CO 2 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.
Our unconstrained ensemble yields a very large spread in the four target variables. δ 13 C DIC ranges from −0.77 ‰ to 1.25 ‰, δ 13 C atm from −0.9 ‰ to 1.32 ‰, CO 2 from −107 to 181 ppm, and CO 2− 3 from −49 to 42 µmol kg −1 . 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 δ 13 C atm values, for instance, is almost uniform within the defined observational range (Fig. 6b), while CO 2 , CO 2− 3 , and δ 13 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 cycleclimate processes as tested within Bern3D. For example, it remains difficult to model a PI-LGM CO 2 of ∼ 90 ppm under the constraint of meeting the other three observational targets. Most ensemble members represent a CO 2 of ∼ 82 ppm (Fig. 6c). Higher CO 2 values are present in the unconstrained ensemble but are generally associated with even more negative CO 2− 3 (not shown), which is outside the defined target range for CO 2− 3 . 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, δ 13 C DIC , is sensitive to this restriction within the constrained ensemble and positive land. The distributions of δ 13 C atm , CO 2 , and CO 2− 3 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 δ 13 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, δ 13 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 δ 13 C DIC . Further, solutions with negative land also yield CO 2 and CO 2− 3 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.

Carbon and δ 13 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 −1 to a mean value over the glacial of 11.25 GtC yr −1 . The oceanic carbon inventory increases at the expense of the atmosphere, yielding a higher CO 2− 3 concentration during the glacial. Turning to the deglaciation, CO 2 increases by 27.8 ppm from LGM to PI, while δ 13 C atm and δ 13 C DIC show little change between LGM and PI (Fig. 7). By applying the standard transient forcings alone, neither the δ 13 C constraints in the atmosphere and ocean nor the CO 2 constraint is met, calling for the consideration of additional processes.

Sensitivity of carbon cycle processes to changes in additional deglacial processes
The mean δ 13 C signature of the ocean-atmosphere system changes in response to potential deglacial processes ( Fig. 7ab). For simplicity, we focus on the change in mean ocean δ 13 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 δ 13 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 δ 13 C DIC . The influence of these mechanisms on δ 13 C DIC cannot be ignored when estimating land based on reconstructed changes in δ 13 C.
Reconstructions of the change in δ 13 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 δ 13 C DIC and δ 13 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 δ 13 C atm , while hardly changing δ 13 C DIC . The strong sensitivity of δ 13 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 δ 13 C independently from mean ocean δ 13 C DIC change (see also Fig. C2b).
Considering CO 2 and CO 2− 3 , the changes in investigated mechanisms significantly influence atmospheric CO 2 and deep Pacific carbonate ion concentration (with the ex-ception of air-sea gas transfer rate). A general relationship becomes apparent. A positive change in CO 2 is always accompanied by a negative change in CO 2− 3 and vice versa ( Fig. 7c-d). Thus, the explanation for the observed deglacial increase in CO 2 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 3 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 7. Sensitivity of (a) δ 13 C DIC , (b) δ 13 C atm , (c) CO 2 , and (d) CO 2− 3 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 2 : Petit et al., 1999;Monnin et al., 2001;Siegenthaler et al., 2005;Lüthi et al., 2008;Köhler et al., 2017; δ 13 C atm : Schmitt et al., 2012; CO 2− 3 : Yu et al., 2013;Qin et al., 2017;Luo et al., 2018; δ 13 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.

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 ‰ −1 for the pair CO 2 and δ 13 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 2 or δ 13 C DIC remains outside the observational range. The four processes are effec-tive in varying δ 13 C DIC but relatively ineffective in varying CO 2 . Processes that have a much larger impact on CO 2 than on δ 13 C DIC are related to the CaCO 3 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 δ 13 C DIC target, while variations in coral reef regrowth, rain ratio, and upper ocean remineralization depth are necessary to meet the CO 2 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 3 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 2 . The associated decrease in δ 13 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 δ 13 C atm − CO 2 resemble the δ 13 C DIC − CO 2 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 δ 13 C atm and low sensitivity for CO 2 , whereas coral reef regrowth, rain ratio, and upper ocean remineralization changes exert a large impact on CO 2 and affect δ 13 C atm little (Fig. 8b). In contrast to the δ 13 C DIC − CO 2 target pair (Fig. 8a), the slopes for the δ 13 C atm − CO 2 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 CO 2− 3 ) gen-erally results in a larger drawdown of atmospheric CO 2 (positive CO 2 ). 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 2 , while at the same time hardly changing deep Pacific CO 2− 3 , making it a potentially important process to fulfill the CO 2 target without shifting CO 2− 3 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 δ 13 C DIC − δ 13 C atm target pair (Fig. 8d). The difficulty in simultaneously fulfilling the atmospheric and oceanic δ 13 C constraints with a negative land is becoming obvious. With negative land, δ 13 C atm and δ 13 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 δ 13 C atm and δ 13 C DIC are both moved towards positive values relative to the standard forc-ings (Fig. 8d). Considering the CO 2 and CO 2− 3 constraints (Fig. 8c) implies that negative land can only be offset by increased glacial organic weathering or rain ratio. However, reaching a CO 2 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).

to changes in deglacial processes
In this section, we describe the spatial distribution of changes in δ 13 C DIC . The reference run with standard LGM forcings (Fig. 9a) shows substantial spatial changes in δ 13 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 δ 13 C DIC in Fig. 9a. In particular, the deepening and strengthening of the AMOC during the deglaciation lead to positive δ 13 C DIC changes in the deep Atlantic and to negative δ 13 C DIC changes in the upper Atlantic, consistent with proxy reconstructions.
The sensitivity simulations may be grouped by processes with small and large changes in the spatial pattern of δ 13 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 atmosphereocean 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 δ 13 C (Fig. 9b). As a result of the carbon removal from the coupled system, changes in bottom water concentrations of CO 2− 3 and O 2 may influence sedimentation fluxes of POC and CaCO 3 locally and lead to small local patterns (see, for instance, negative anomalies in δ 13 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 δ 13 C DIC (Figs. 9h and 7a). Increased rates of coral reef regrowth have little impact on δ 13 C DIC as the isotopic signature is similar for DIC and deposited CaCO 3 (Fig. 9f). The slight increase in δ 13 C DIC in response to enhanced coral reef regrowth has been attributed to changes in fractionation during marine photosynthesis (Freeman and Hayes, 1992;. Changes in the other four mechanisms affect the spatial pattern of δ 13 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) andMenviel et al. (2015), a more poorly ventilated ocean due to reduced wind stress leads to an increase in δ 13 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 13 C-depleted carbon from the ocean increases δ 13 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 3 cycling spatially, overlay each other, yielding the δ 13 C DIC distribution seen in Fig. 9C. Ocean-sediment interactions also play an important role for δ 13 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 δ 13 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 3 increase in response to reduced oxygen concentrations from a lower gas transfer rate and changes in the CO 2− 3 concentration, as does the δ 13 C signature of the flux, leading to a removal of 13 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 δ 13 C DIC still prevails (Fig. 9d).
As discussed in detail in Tschumi et al. (2011), responses in δ 13 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 3 and POC. This negative perturbation in δ 13 C is not removed entirely over the deglacial and Holocene due to the long ocean-sediment response timescale, leaving the generally negative δ 13 C DIC as seen in Fig. 9e. 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 δ 13 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).
in δ 13 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 δ 13 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 δ 13 C signatures, which are then reflected in POC and CaCO 3 export fluxes and coral reef regrowth, and altered sedimentation fluxes as a result of altered bottom water CO 2− 3 and O 2 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 oceansediment interactions shapes the patterns in δ 13 C DIC .

Discussion
We show that carbon exchange with ocean sediments and the lithosphere biases earlier δ 13 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 pulserelease 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 δ 13 C signature significantly. Their influence cannot be ignored when analyzing the budget of δ 13 C on deglacial timescales. The importance of ocean sediment and lithosphere fluxes in regulating past atmospheric CO 2 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., , 2018 reveal that these fluxes are more dynamic than previously assumed.

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 2 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 3 , 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 δ 13 C of DIC agree with corresponding reconstructions (see Fig. 9 in . 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 2 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 δ 13 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 δ 13 C (Broecker, 1982;Zimov et al., 2006Zimov et al., , 2009Zech et al., 2011;Kemppinen et al., 2018). The deglacial change in CO 2 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 −1 at the LGM to around 17 GtC kyr −1 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 866 A. Jeltsch-Thömmes et al.: Low terrestrial carbon storage at the Last Glacial Maximum 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.

Contribution of deglacial mechanisms to the CO 2 rise
Turning to the role of different mechanisms for the deglacial increase in atmospheric CO 2 , we find that ocean circulation changes only partly explain the deglacial CO 2 rise. Ocean circulation changes are not able to simultaneously explain glacial-interglacial changes in δ 13 C DIC and CO 2 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 2 , while at the same time hardly changing deep Pacific CO 2− 3 . 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 . This renders this class of processes potentially important to explain part of the deglacial CO 2 increase, as reconstructions suggest small LGM to PI changes in deep Pacific CO 2− 3 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 alkalinitybased mechanisms for the deglacial CO 2 increase (see Fig. C2d). The constrained ensemble yields a large extra burial of CaCO 3 as well as an increase in the rain ratio of CaCO 3 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 3 burial flux (below 200 m) increased by about 0.2 GtC yr −1 from the LGM to the Holocene.

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(Zimov et al., , 2009Zech 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(Zimov et al., , 2009Zech 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 δ 13 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 . Ice core CO 2 and δ 13 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.

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 δ 13 C signature of atmospheric CO 2 and of dissolved inorganic carbon in the ocean; and the deglacial change in atmospheric CO 2 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 13 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 13 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 13 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 2 . 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 2 . The results demonstrate that ocean sediments and the weathering-burial cycle are an integral part of the Earth system playing a fundamental role on glacialinterglacial timescales. 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 2 and 14 CO 2 are implemented according to OCMIP-2 protocols (Najjar and Orr, 1999;Orr et al., 1999) with updates for the 14 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 14 C, the global mean air-sea transfer is reduced by 19 % compared to OCMIP-2 in order to match observation-based radiocarbon estimates . 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 13 CO 2 .
Biogeochemical cycling in the model is detailed in Tschumi et al. (2011) and with further documentation of results in follow-up studies (e.g., Menviel et al., , 2015Roth and Joos, 2012;Roth et al., 2014;Battaglia et al., 2016;Battaglia and Joos, 2018b, a). Dissolved inorganic carbon and semilabile organic carbon (DIC, DOC), the corresponding isotopic forms (DI 13 C, DO 13 C, DI 14 C, DO 14 C), and alkalinity (Alk), phosphate (PO 4 ), oxygen (O 2 ), 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 4 , Alk, and O 2 and the corresponding elemental ratios in DOC and POC are linked by fixed Redfield ratios (P : Alk : C : O 2 : Alk = 1 : 17 : 117 : −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 where z 0 refers to the reference depth of 75 m (see also Roth et al., 2014). Export of calcium carbonate (CaCO 3 ) and of opal is computed from new production and availability of dissolved silica in the euphotic zone. The export rain ratio of CaCO 3 to POC is set to 0.083 in the absence of dissolved silica, while the production ratio of CaCO 3 to POC is reduced at the favor of opal production in regions with abundant availability of silicic acid. Particulate CaCO 3 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 3 , opal, POC, clay) and seven tracers in the pore water (DIC, DI 13 C, DI 14 C, alkalinity, phosphate, oxygen, and silicic acid) are modeled. The dissolution of CaCO 3 depends on the pore-water CO 2− 3 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 2 is not consumed below a threshold, somewhat reflecting the process of denitrification without modeling NO − 3 . The respective reaction rate parameters for CaCO 3 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 4 , Alk, DIC, DI 13 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 2 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 13 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 13 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. 13 C stocks and flows are modeled in the atmosphereocean-sediment-land biosphere system. 13 C fluxes to the lithosphere associated with the burial of POC and CaCO 3 and weathering input are explicitly simulated. 13 C fractionation is considered for air-sea gas transfer, carbonate chemistry, the formation of CaCO 3 , 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 13 CO 2 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 3 formation is small and computed following Mook (1986). This results in an isotopically heavy signature of around 2.9 ‰ for CaCO 3 . Fractionation during photosynthesis in the ocean and thus between dissolved CO 2 ([CO 2 ]) and POC and DOC is calculated according to Freeman and Hayes (1992). The fractionation increases logarithmically with [CO 2 ] in surface water and results in an isotopically light signature for POC and DOC of around −20 ‰. A lowering of [CO 2 ], or correspondingly of pCO 2 , by about 10 % yields a change in the isotopic signature of POC by +0.55 ‰. The burial of POC and CaCO 3 results in an isotopic signature of the burial flux of around −9 ‰, intermediate between the isotopically light POC and the isotopically heavy CaCO 3 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 2 is set to 278 ppm and δ 13 C of CO 2 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 2 concentrations following the parameterizations given in Colbourn et al. (2013), accounting for changes in CaCO 3 and CaSiO 3 weathering on land.
Simulations are started from the end of the PI spin-up, and atmospheric CO 2 and δ 13 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 δ 13 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 2 and δ 13 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: 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.

870
A. Jeltsch-Thömmes et al.: Low terrestrial carbon storage at the Last Glacial Maximum

Appendix B: Sediment tuning
The aim of the sediment tuning was to improve the representation of export, deposition, and dissolution fluxes of POC, CaCO 3 , 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 3 in the reactive uppermost 10 cm of the sediment represented in the sediment module. Prior to tuning, the global CaCO 3 content in the model was ∼ 400 GtC. Using observational constraints on modern global bulk mass accumulation rates (MARs) , 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 3 stock of ∼ 1000 GtC. The second concern was the opal burial that used to be ∼ 3 Tmol Si yr −1 but with new observational data suggesting a burial flux of ∼ 6.5 Tmol Si yr −1 (Tréguer and De La Rocha, 2013). Further, all other sediment-related fluxes of POC, CaCO 3 , and opal were tuned with regard to observational constraints (see Table 1). 2900 m 5066 m Diffusion coefficient in the sediment pore water 8 × 10 −6 cm 2 s −1 5.5 × 10 −6 cm 2 s −1 Reactivity for calcite in the sediment 1000 L mol −1 yr −1 800 L mol −1 yr −1 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 100member 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 Zweng et al., 2013;Garcia et al., 2013a, b) datasets were used. For export, deposition, and burial of POC, CaCO 3 , and opal values from Battaglia et al. (2016), Tréguer and De La Rocha (2013), Sarmiento and Gruber (2006), Milliman andDroxler (1996), andFeely et al. (2004) were used. For the sediment data compilations by Cartapanis et al. (2016Cartapanis et al. ( , 2018 were used to compare the model performance. The distribution of CaCO 3 in the sediments after the tuning compared to observational data is shown in Fig. B1. Export, deposition, and burial fluxes for POC, CaCO 3 , 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).  (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.