Late Pleistocene Carbon Cycle Revisited by Considering Solid Earth Processes

The importance of volcanic CO2 release, continental weathering, and coral reef growth on the global carbon cycle has been highlighted by several different studies. Based on these independent approaches, we here revisit the last 800 kyr with the box model BICYCLE, which has been extended to be able to address these solid Earth contributions to the carbon cycle. We show that the volcanic outgassing of CO2 as a function of sea level change from mid‐ocean ridges and hot spot island volcanoes cannot be the generic process that leads during phases of falling obliquity to a sea level‐CO2 decoupling as has been suggested before. The combined contribution from continental and marine volcanism, if both lagging sea level change by 4 kyr, might have added up to 13 ppm to the glacial/interglacial CO2 rise. Coral reef growth as suggested by an independent model is during glacial terminations about an order of magnitude too high to be reconciled with meaningful carbon cycle dynamics. Global riverine input of bicarbonate caused by silicate and carbonate weathering is suggested to have been stable over Termination I. However, if weathering fluxes are changed by up to 50% in sensitivity experiments, the corresponding bicarbonate input might contribute less than 20 ppm to the deglacial atmospheric CO2 rise. The overall agreement of results with the new process‐based sediment module and the previously applied time‐delayed response function to mimic carbonate compensation gives confidence in the results obtained in previous applications of the BICYCLE model without solid Earth processes.


Introduction
Atmospheric CO 2 is a global integrating variable that captures changes in various different parts of the global carbon cycle. On glacial/interglacial (G/IG) timescales of the late Pleistocene, the period for which reliable reconstruction of CO 2 from ice cores exist (Bereiter et al., 2015), most hypotheses trying to explain parts of the observed CO 2 changes deal with marine processes (e.g., Ganopolski & Brovkin, 2017;Khatiwala et al., 2019), since small changes in the ocean that contains nearly 2 orders of magnitude more carbon than the atmosphere are easily able to explain rather larger changes in the atmospheric CO 2 content. This oceanic-centered view is not wrong, and variability in ocean circulation and the marine carbon pumps cannot be ignored. However, results are still highly scenario-and model-dependent (e.g., see Gottschalk et al., 2019), simply for the fact that the atmosphere is-simply spoken-the gas phase of the surface ocean, and atmospheric CO 2 is most directly determined by surface ocean alkalinity, dissolved inorganic carbon (DIC), and gas exchange processes.
During the last decade or so, evidence from both changes in terrestrial carbon storage and of fluxes between the solid Earth and the atmosphere-ocean subsystem have been suggested to play a larger role for late Pleistocene carbon cycle dynamics. The challenge of reconstructing especially the build-up and release of carbon from terrestrial reservoirs with long turnover times, such as peatlands and permafrost, has been facilitated by some data-based suggestion (e.g., Ciais et al., 2012;Meyer et al., 2019;Winterfeld et al., 2018) but is still far from been completely understood. We therefore do not investigate in that direction any further but note that others have done so (e.g., Ganopolski & Brovkin, 2017).
The focus of this study will be on suggested changes in solid Earth carbon fluxes that might be of interest on G/IG timescales. In detail, we have a closer look on volcanic CO 2 outgassing and weathering. A precondition for such studies is the implementation of a process-based model of ocean-sediment exchange of fluxes of CaCO 3 in any carbon cycle model applied. Various carbon cycle models of different complexity have implemented such ocean-sediment interactions (e.g., Brovkin et al., 2007;Heinze et al., 1999; Furthermore, coral reef production over the Quaternary has been estimated from a modeling approach (Husson et al., 2018), leading to annual CaCO 3 fluxes in shallow waters, whose consequence for the global carbon cycle as a whole needs to be tested. In detail, Husson et al. (2018) find that coral reefs produced up to 19 Tmol/yr of CaCO 3 during Termination I and the Early Holocene, more than 4 times as much as the previous estimate by Vecsei and Berger (2004). Previous model applications of shallow water accumulation of CaCO 3 in coral reefs are, for example, Ridgwell et al. (2003) or Menviel and Joos (2012). Alternatively, a simpler approach of coral reef growth as function of sea level change might be taken into consideration (Munhoven & François, 1996).
For this study we implement the independently suggested changes of those processes (volcanic CO 2 outgassing, continental weathering, and coral reef growth) in a revised version of the carbon cycle box model BICYCLE ) that has the advantage to run very fast and that also simulates changes in 13 C and 14 C and discuss resulting consequences. BICYCLE has furthermore been extended by an ocean-sediment exchange module (Munhoven & François, 1996;Munhoven, 1997), consisting of a mixed sedimentary zone on top of a variable number of "historical" layers of sediment that can now-instead of a so far used simplified feedback mechanism of carbonate compensation-simulate the flux of CaCO 3 from the ocean to the sediments-and, if the saturation boundary conditions change, also the reverse flux following the dissolution CaCO 3 in the sediments.

Brief Summary of the General Model Setup of BICYCLE
The core of BICYCLE (the Box model of the Isotopic Carbon cYCLE) is a 10-box ocean (O) and a 7-box terrestrial biosphere (B) together with a 1-box atmosphere (A), in which the concentration of carbon (as DIC in the ocean, as pCO 2 in the atmosphere, as organic carbon in the biosphere) and both of the isotopes 13 C and Δ 14 C are traced (Köhler et al., 2005). Furthermore, in the ocean alkalinity, PO 3− 4 as macro-nutrient and O 2 is represented. From the two variables of the marine carbonate system (DIC and alkalinity) all other variables (CO 2 , HCO − 3 , CO 2− 3 and pH) are calculated (Zeebe & Wolf-Gladrow, 2001) with updates of the dissociation constants pK 1 and pK 2 (Prieto & Millero, 2002). In detail, all parts of the marine carbonate system are calculated as function of temperature, salinity and pressure. The distribution of temperature, salinity and the various tracers have been initialized for pre-industrial conditions from the World Ocean Atlas, while changes are calculated following the applied forcing described below. The net gas exchange between the atmosphere and the surface ocean reservoirs is driven by the difference between the partial pressure of CO 2 in the atmospheric box and that in the surface water, derived from the concentration of CO 2(aq) using Henry's law, where the solubility of CO 2 is calculated as a function of temperature and salinity. In addition, we use a piston velocity of 2.5 and 7.5 × 10 5 m s −1 for low and high latitude reservoirs, respectively. The average gas exchange coefficient (0.051 mol m −2 yr −1 μatm −1 ) is of the same order as in carbon isotope studies (Heimann & Maier-Reimer, 1996) and leads to pre-industrial gross fluxes of CO 2 between ocean surface and atmosphere and vice versa of about 60 Pg C/yr. The pre-industrial marine export production of organic carbon at 100 m water depth is set to 10 Pg C/yr with a fixed molar rain ratio of organic C/CaCO 3 of 10:1. Export production depends via the applied elemental (Redfield) ratios of C:N:P:O 2 = 123:17:1:−165 (Körtzinger et al., 2001) on available macro-nutrients (here PO 3− 4 >) and is realized first in the equatorial regions, and filled up to the prescribed numbers in polar regions thereafter. These fluxes of the marine biology have been taken from full ocean general circulation models (OGCMs) that have been validated with field data (e.g., Sarmiento et al., 2002;Schlitzer, 2002). From the marine biological production in the surface layers, 28% of the organic material and 80% of CaCO 3 reach the deep ocean below 1,000 m depth. For deep ocean O 2 concentration <4 μmol/kg 3 the remineralisation of organic matter is assumed to follow the denitrification pathway and thus does not consume any molecular oxygen, in line with the Ocean Carbon-Cycle Model Intercomparison Project (OCMIP) 2 protocol. In the previous versions of BICYCLE the fluxes of CO 2− 3 between ocean and sediment have been simplified as carbonate compensation feedback to any changes in oceanic CO 2− 3 , as described in Köhler and Fischer (2006). Prescribed ocean circulation fluxes have been revised before in order to meet data constraints given by 13 C . The seven-box terrestrial biosphere (Köhler & Fischer, 2004) distinguishes between the C 3 and C 4 photosynthetic pathways of grasses, considers non-woody and woody parts of trees, and contains three soil carbon boxes with turnover times ranging from 5 to 1,000 years. Terrestrial carbon storage depends via CO 2 fertilization directly on atmospheric CO 2 , and via changing respiration rates on global temperature. A plot of the different boxes within this AOB model configuration is shown in Figure 1b. The AOB configuration of BICYCLE does not consider volcanic CO 2 emissions, weathering fluxes, or coral reef growth.

Time-Dependent Forcing of the Model Without Solid Earth Processes
The previous applications of BICYCLE to climate changes of the Late Pleistocene have been restricted to the last 740 kyr due to data availability from ice cores (e.g., Köhler & Fischer, 2006). Here, we apply the model for the last 800 kyr, which represents the time window of nowadays available ice core data. This extension in application time leads to some adjustments in the forcing data sets: The EPICA Dome C (EDC) D ice data set (Figure 2c), which has been taken to determine temperature changes in the Southern Ocean, is now available in its full length (EPICA-community-members, 2004;Jouzel et al., 2007). Here, D ice is corrected for sea level changes via 18 O sw of sea water (sw) following D corr = D ice − 8 ⋅ 18 O sw ⋅ (1 + D ice ∕1,000)∕(1 + 8 ⋅ 18 O sw ∕1,000)  using 18 O sw from Bintanja et al. (2005).
However, the EDC iron flux (Wolff et al., 2006), which has been taken by BICYCLE to force enhanced marine biological productivity in the Southern Ocean in glacial times, is still only available for the last 740 kyr. We therefore substitute this forcing data set by dust fluxes (Lambert et al., 2008) and similarly assume that the threshold value of this variable at 18 kyr BP determines when the Southern Ocean marine biology switches from an iron-unlimited to an iron-limited regime. In the iron-unlimited case potentially all available surface PO 3− 4 in the Southern Ocean might be used for export production, rising the global export of organic matter at 100 m water depth to 13 Pg C/yr, while in the iron-limited case this global export production is limited to 10 Pg C/yr ( Figure 2d).
All ice core data have been transferred to the most recent age model AICC2012 Veres et al., 2013). The full used forcing data sets have been plotted in Figure 2 in Köhler and Fischer (2006). We here show the updated and extended ice core data from EDC together with most other important forcing records and the consequential imposed temporal changes on the boundary condition in Figure 2 including sea level, mean ocean salinity, North Atlantic Deep Water (NADW) formation strength; Southern Ocean vertical mixing, global integrated export production of organic matter at 100 m water depth ocean sea surface, deep and mean ocean temperatures together with global integrated sea ice area. . V: outgassing of CO 2 from volcanoes on land potentially and temporally overlain by land ice and from hot spot island volcanoes (and mid-ocean ridges, not shown) influenced by changing sea level; C: shallow water carbonate deposition due to coral reef growth; Si-W: silicate weathering and Ca-W: carbonate weathering with different sources of C, but both delivering HCO − 3 ions into the ocean; P: PO 3− 4 riverine input and sedimentary burial; S: CaCO 3 sedimentation and dissolution. A2B: atmosphere-biosphere exchange of CO 2 ; A2O: atmosphere-ocean exchange of CO 2 . The cyan-colored broken circles mimic the two overturning cells in the Atlantic and Indo-Pacific Ocean. The previous model version has been restricted to the fluxes A2O, A2B, and some simplified representation of carbonate compensation between ocean and sediment. (b) Details of subsystem AOB  with the prescribed pre-industrial water fluxes in the ocean (blue, numbers in Sv). Black arrows denote carbon fluxes. The terrestrial biosphere contains: C4: C 4 ground vegetation; C3: C 3 ground vegetation; NW: non-woody parts of trees; W: woody parts of trees; D: detritus or above-ground litter; FS: fast decomposing soil or below-ground litter; SS: slow decomposing soil. (c) A conceptual sketch of the new sediment module, the figure is taken from Munhoven (1997). Z j is a running index with j ∈ [0, 80] addressing the water depth (layers of 100 m each) at which sedimentary processes in each water basin are considered.
In detail, the interglacial ocean circulation depicted by blue arrows in Figure 1b has been prescribed from WOCE data (Ganachaud & Wunsch, 2000) showing a NADW formation strength of 16 Sv. Based on benthic 18 O in the sediment core ODP980 in the North Atlantic and the detailed analysis of the carbon cycle dynamics during Termination I (Köhler et al., 2005) a threshold in 18 O is defined for the switch between interglacial and glacial NADW formation strength (Figure 2b), the latter being 10 Sv based on various modeling studies (e.g., Meissner et al., 2003). Millennial-scale variability in Atlantic meridional overturning circulation related to the bipolar seesaw (e.g., Stocker & Johnsen, 2003;Zhang et al., 2014) are ignored here. The vertical water mass exchange flux from the surface to the abyss in the Southern Ocean was coupled linearly to Southern Ocean temperature changes (which are themselves prescribed from the EPICA Dome C D  (Flower et al., 2000;McManus et al., 1999;Wright & Flower, 2002), dotted line marks the threshold for switching between both states. (c) EDC ice core δD (EPICA-community-members, 2004;Jouzel et al., 2007) corrected for δ 18 O sw , from which Southern Ocean sea surface temperatures (SST) and vertical mixing (SO-x: SO surface-to-deep ocean flux, right y axis) is calculated. (d) Marine biology in the Southern Ocean (SO) is either Fe-limited or Fe-unlimited following dust fluxes in the EDC ice core (Lambert et al., 2008). The dotted line marks the threshold for switching between both states, leading to global integrated export production of organic matter at 100 m water depth (right y axis). (e) Different ocean temperatures averaged within the model. SST is calculated from all ocean surface boxes, deep ocean temperature from the boxes with water depths below 1,000 m. The sea ice free SST is relevant for air-sea gas exchange. Right y axis: global integrated sea ice area. In Subfigures (b)-(d) original data (thin lines) and 3-kyr running mean (bold lines) are shown. BICYCLE is forced by the running-mean data.
record) with 29 Sv during preindustrial times and 9 Sv during the LGM (Figure 2c). The reversed upwelling flux in the Southern Ocean changed identically to the downwelling flux. Apart from the overturning cell in the Atlantic that mimics the Atlantic overturning circulation (and which is directly related to NADW formation strength), the Southern Ocean vertical mixing fluxes, a closure of the Bering Strait in the glacial setting and some necessary readjustments to allow water mass conservation all other water fluxes are kept at their values depicted in Figure 1b. The glacial circulation regime that actually was applied during most of the times during the last 800 kyr is plotted in Figure 1 of Köhler et al. (2010). Sea ice in the high latitudinal areas was prescribed for interglacial values from observations (Cavalieri et al., 1997) to 20 × 10 12 m 2 globally. Past changes in fractional sea ice cover in the high latitudinal surface boxes are based on LGM reconstructions (Gersonde et al., 2005;Sarnthein et al., 2003) leading to a global LGM sea ice cover of 36 × 10 12 m 2 ( Figure 2e). Based on these two tie points, all other sea ice areas are calculated linearly from changes in local SST. Further details on model forcing are found in Köhler and Fischer (2006).
In previous applications, we had considered a series of forcings that could be switched on or off in order to analyze their impacts on the global carbon cycle changes. These included ( previous applications of the model the contribution of individual processes to changes in atmospheric CO 2 will also be briefly discussed at the end of this study.

Model Extension Toward Solid Earth Processes
The model extension toward the version BICYCLE-SE used here is outlined in Figure 1a. The main improvement is the implementation of a sediment module that captures early diagenesis in the sedimentary mixed layer, under which there is a pile of historical layers making up a synthetic sediment core (Figure 1c). Implementation and realization of the sedimentary processes directly follows Munhoven and François (1996) and Munhoven (1997). In each mixed and historical layer in the sediment the amount of CaCO 3 , of other material (clay) and a time marker recording the average year of sedimentation is traced. This time marker is then used to calculate an offset in Δ 14 C from deep ocean values once the mixed-layer is subject to chemical erosion and the most recently deposited historical layer gets mixed back into the sedimentary mixed-layer again. The carbon isotope 13 C in the sedimentary mixed layer are only followed in detail in aggregated boxes (one box for each of the three ocean basins), since this setup had already been implemented in the model and used in previous applications within the AOB configuration.

Details on the Sediment Module
In each basin the carbonate chemistry in the water column is calculated for different water depths in steps of z = 100 m. We are thus able to consider the pressure-dependency of the ocean-sediment interaction, even if the concentrations of DIC and alkalinity stay constant below 1,000 m water depth due to our box model setup. For each of these different water depths a basin-wide average mean sediment column is simulated consisting of a sedimentary mixed layer of 8 cm on top of various historical layers, each 8 mm thick (Figure 1c). The areal extend of such a sediment column at a specific water depth z in the different oceanic reservoirs is calculated following the geometry based upon realistic bathymetric profiles. The ocean-sediment exchange happens only via this sedimentary mixed layer. All CaCO 3 exported from the surface ocean is together with other material (e.g., clay) accumulated in the sedimentary mixed layer. The amount of other material is calculated based on a fixed mass ratio of 0.1 (other/CaCO 3 ) using a density of 2,900 and of 3,000 kg/m 3 for CaCO 3 and others, respectively. Furthermore, we assume a porosity of the sediments of 0.7. The total CaCO 3 export is partitioned into an aragonite and a calcite fraction with 65% being aragonite. These parameter values have been taken unadjusted from Munhoven (1997). The fraction of aragonite at CaCO 3 is within the range of 33-89% of a most recent estimate (Buitenhuis et al., 2019) but higher than the 35% found in Gangstø et al. (2011). Sensitivity tests have shown a value of the order of 35% would lead to a drop in deep ocean CO 2− 3 by 10-15 μmol/kg and to a rise in atmospheric CO 2 by 20 ppm. Both changes would bring simulated carbon cycle values away from reconstructions and such alternative values for the aragonite share in CaCO 3 export are therefore neglected. Aragonite that reaches sedimentary reservoirs located above the current Aragonite Saturation Horizon (ASH) is transferred to the sedimentary reservoirs and "converted" to calcite; aragonite that falls deeper than the ASH is instantaneously dissolved and not even transferred to the sediment boxes. We thus adopt a semi-reflective bottom boundary condition for aragonite.
If the sedimentary mixed layer increases in thickness by more than 60% of a historical layer of 8 mm, a new historical layer is established, reducing subsequently the amount of material in the sedimentary mixed layer accordingly. In water depths deeper than the calcite saturation horizon (CSH), the carbonate chemistry allows the partial dissolution of CaCO 3 in the sediments following early diagenesis which is proportional to the carbonate ion undersaturation to the power of (n + 1)/2 with n = 4.5 being the order of the dissolution law (Keir, 1980). More details on the sediment module are found in Munhoven (1997), but the finally implemented equation for early diagenesis is contained in Appendix A1. The sedimentary mixed layer might thus become thinner. Once its thickness drops by more than 60% of an historical sediment layer the layer directly underneath the mixed zone is included into the mixed layer. In effect, the sedimentary mixed layer varies in thickness between 7.52 and 8.48 cm (±60% of an historical layer of 8 mm around its standard depth of 8 cm).
For the initialization of the sediment variables (CaCO 3 content, time since sedimentation) the model is run in its solid Earth setup and the derived values of the sedimentary variables after 800 kyr are stored and used as initial values for all further simulation scenarios.

Implementing Solid Earth Processes
Once a sediment module is implemented in a carbon cycle model one has to also to consider other solid Earth fluxes, otherwise the simulated system quickly drifts away due to an imbalance of carbon sources and sinks. Thus, with the aim to avoid long-term drifts in the global means of alkalinity, PO 3− 4 , DIC, CO 2 , Figure 3. Results from Scenario SE with focus in ocean-sediment interaction (solid lines) compared to some output from Scenario AOB without sediment module and without solid Earth contributions (long broken lines). From the sediment module the carbonate fraction (shaded area covers depths from 10% [lower bound] to 90% [upper bound] with 50% as short broken line in between) and the calcite saturation horizon (CSH) in the water column as a function of water depth are plotted for the equatorial (a) Atlantic and (b) Indo-Pacific. Furthermore, for both ocean basins the deep ocean CO 2− 3 concentration and its reconstructions for the last 160 kyr (gray from Yu et al., 2013) and in the Pacific for the last 250 kyr (yellow from Qin et al., 2017) are plotted. Deep Atlantic data from core VM28-122 (12 • N, 79 • W, 1,800 m sill depth), deep Pacific data from PC61 (1 • S, 140 • W, 4,276 m water depth), and the mean from two species-specific reconstructions from core WP7 (yellow, 4 • S, 156 • E, 1,800 m water depth). (c) Total alkalinity (TA), sum of carbon in AOB subsystem, in the ocean (O). (d) Global integrated carbon in the sedimentary mixed layer (M) in Scenario SE. (e) Atmospheric CO 2 in SE and AOB against ice core data (gray circles from Bereiter et al., 2015). r 2 are calculated from linear regression against the data. (right axis) Atlantic-Pacific difference in deep ocean CO 2− 3 for the same simulations and data shown in Subfigures (a) and (b). The output from the sediment module has been generated every 1 kyr, while all other data have been written out every 100 years. The vertical broken line marks the LGM (20 kyr BP). and 13 C we derive a new standard run (Scenario SE). It consists of a constant background level of volcanic CO 2 outgassing (∼9 ⋅ 10 12 mol C/yr) on which some variability that depends on changes in sea level is added (volcanic Scenario V3), Furthermore, Scenario SE contains constant weathering fluxes of 12 ⋅ 10 12 mol/yr HCO − 3 riverine input for silicate and carbonate weathering each (weathering Scenario W0), and a simplistic coral reef growth as function of sea level rise (coral reef Scenario C1).
To avoid long-term drifts in 13 C a sedimentary burial of a fraction of the organic carbon that reaches the abyss-and the related nutrient PO 3− 4 -is now implemented. The fixed fraction of ∼0.6% of the export production of organic carbon that reaches the deep ocean boxes is permanently buried in the sediment, representing 35 ⋅ 10 12 g C/yr. The exact amount has been determined by avoiding long-term trends in 13 C, since the loss of organic C that is with approximately −15‰ depleted in 13 C is counterbalancing other new fluxes: carbonate weathering introduces C with a 13 C of +2‰ (an average of carbonate rocks) and volcanic CO 2 Table 1 Overview of Simulation Experiment Scenarios Acronym Description

AOB10
Previous CTRL run , 740 kyr long, AOB with carbonate compensation as simple response function AOB 800 kyr long with updated forcing, all else as AOB10
with on average a 13 C of −3.8‰ (de Leeuw et al., 2007;Huybers & Langmuir, 2009;Sano & Williams, 1996;Shaw et al., 2003). Furthermore, carbonate lost to the ocean due to coral reef growth leaves the ocean with a 13 C of −1.36‰ which is the mean from 41 pre-Suess effect samples in the review of Linsley et al. (2019).
With the chosen Redfield ratio in organic matter of C:P = 123:1 (Körtzinger et al., 2001) the sedimentary loss of PO 3− 4 is balanced by a constant riverine input of 24 ⋅ 10 9 mol P/yr. A riverine input of organic matter is not envisaged, since Regnier et al. (2013) showed that any carbon flux from land to the open ocean is to a large extend generalized into DIC or lost to either sediment or atmosphere during the passage through inland waters, estuaries and coastal ocean. Note that this organic carbon sink in the sediment is only implemented to satisfy 13 C but not as a process in the sediment module that might influence the dissolution of CaCO 3 . In the discussion section we will expand on further motivations for this assumption and how this might influence our simulation results. The consequences of this solid Earth setup in comparison to the AOB-only version for the ocean-sediment exchange processes are plotted in Figure 3, which will be discussed in detail in the results section.
Equipping BICYCLE with a process-based sediment module enables the revised model version BICYCLE-SE to address questions related to solid Earth carbon fluxes (volcanic outgassing of CO 2 , weathering of silicate or carbonate rocks on land, and coral reef growth as a shallow-water carbonate sink). Implemented processes and our motivation for chosen setups are given in detail in the following subsections. All simulation scenarios are summarized in Table 1, and the numerical realization of some of the scenarios is compiled in the Appendix B.

CO 2 Release From Volcanism
Volcanic outgassing of CO 2 is known to be one of the main long-term drivers of the global carbon cycle and climate over multi-million years timescale (e.g., Berner, 1990). The volcanic carbon emissions into the atmosphere need on average to be balanced by the removal from continental weathering. In order to avoid drifts in the oceanic DIC and alkalinity budgets the long-term average (say over 10 4 -10 5 years) volcanic CO 2 emission rate and the long-term average CO 2 consumption rate by silicate weathering must be in a ratio of 1:2 (Munhoven & François, 1996). Since the burial of organic carbon in the sediment was also considered here (this was ignored in the analytical calculation of Munhoven & François, 1996), this ratio needs to be slightly adapted.
In addition to such long-term effects it has now become evident that the release of ice sheet mass load might increase subaerial volcanic activity of continental volcanoes and CO 2 outgassing (Huybers &   . Impact of volcanic CO 2 outgassing on atmospheric (a, b) CO 2 . (c) Different volcanic CO 2 emission scenarios (V0, V1, V2, V3) and sea level, on which they partly depend. V0: volcanic CO 2 outgassing constant; V1: additional contribution due to the release of land ice masses on the continents (Huybers & Langmuir, 2009); V2: additional contribution due to the release of sea level on mid-ocean ridges and hot spot island volcanoes (Hasenclever et al., 2017); V3: merging effects of V1 and V2. Gray points (±1σ) in (a) refer to ice core reconstructions (Bereiter et al., 2015). In Subfigure (c) the fluxes are plotted as 1-kyr running mean. Gray line in (c) is the imposed sea level (Bintanja & van de Wal, 2008). Langmuir, 2009). This additional CO 2 outgassing is of relevance on G/IG timescales and lags the ice sheet decrease by ∼4 kyr (Kutterolf et al., 2013), and should also work the other way round (an increase in ice load reduces the volcanic CO 2 fluxes). The magnitude of this additional CO 2 source is highly uncertain but has been estimated to provide on average 8-9 ⋅ 10 12 mol C/yr during Termination I (Huybers & Langmuir, 2009;Roth & Joos, 2012). We implement such a flux as a function of sea level change, which serves as a proxy for the land ice volume change (Scenario V1) that leads (time-delayed by 4 kyr) to a CO 2 flux of at maximum 480 ⋅ 10 12 mol C per m of sea level rise (caused by land ice shrinking), comparable with the deduced full amplitude in the initial study of Huybers and Langmuir (2009).
Furthermore, Hasenclever et al. (2017) have shown that a reduction of water mass in the oceans due to a fall in sea level might also has consequences for CO 2 outgassing from hydrothermal mid-ocean ridges and hot spot island volcanoes. A similar idea has been proposed by Huybers and Langmuir (2017) but focusing only on CO 2 outgassing from mid-ocean ridges. This sea level effect proposed by Hasenclever et al. (2017) works in principle the opposite way than that of Huybers and Langmuir (2009). We therefore adopt here (Scenario V2) a response of the marine volcanic system with an assumed magnitude of 160 ⋅ 10 12 mol C per m of sea level fall as calculated by solid Earth models in Hasenclever et al. (2017). This is equivalent to one third (and reverse) of the amplitude related to land ice mass load considered in Scenario V1. We furthermore simplify the considered processes by assuming that all emissions are subaerial here, while Hasenclever et al. (2017) have shown that nearly half of them might indeed be submarine. However, first carbon cycle experiments in the same study have shown that the differences of subaerial and submarine emissions for atmospheric CO 2 are negligible on the long run; it is only at peak height that they differ by about 1 ppm.
While Hasenclever et al. (2017) focused on the effect of sea level fall, we again assume that this process works the other way around for reversed forcing, leading to a reduction in volcanic emissions during sea level rise, and thus working in opposite direction of the observed CO 2 concentration obtained from ice core data. In the final Scenario V3, which is also chosen for the standard run SE, both processes related to the land ice (V1) and to sea level (V2) are in operation together ( Figure 4c). As response time between sea level change and hydrothermal activity we here assume for simplicity again 4-kyr-the same time delay as for land ice volume change and continental volcanism. This implies that in Scenario V3 both fluxes change simultaneously, but opposite to each other, as a function of sea level change, and thus add up to a minimal contribution of volcanism to atmospheric CO 2 : The sea level-related component simply reduces the land ice-related contribution. We discuss below how CO 2 might change if the response times of both processes are different. However, we acknowledge that our knowledge on this Paleoceanography and Paleoclimatogy 10.1029/2020PA004020 response time between sea level change and hydrothermal activity is limited. Data from iron flux variability suggest an instantaneous response, at least of iron release at mid-ocean ridges and sea level change (Middleton et al., 2016). Hasenclever et al. (2017) summarized that the time lags between change in sea level and a hydrothermal response might be 5-15 kyr, while Huybers and Langmuir (2017) suggested (largely based on Burley & Katz, 2015) a time delay of 50 kyr. The frequency analysis of the crustal thickness from the East Pacific Rise has a maximum cross-correlation with sea level when lagging the latter by ∼45 kyr (Boulahanis et al., 2020). Another model-based study, however, shows that peaks in volcanic emissions from mid-ocean ridges lag peaks in sea level by less than 20 kyr (Cerpa et al., 2019). Lund and Asimow (2011) find response times between sea level change and mid-ocean ridge hydrothermal activity between less than 1 and 15 kyr.
Our constant flux of background volcanic activity with a CO 2 emission of about 9 ⋅ 10 12 mol C/yr falls in the upper part of the range of fluxes compiled (1-12 ⋅ 10 12 mol C/yr) in the review of Burton et al. (2013), in which the upper end refers to the best guess of the authors. A recent review (Fischer & Aiuppa, 2020) derives a total sum of volcanic subaerial CO 2 emissions that is with 3-4 ⋅ 10 12 mol C/yr smaller than our value, but this review also denotes that only Burton et al. (2013) covers the most comprehensive attempt of volcanic CO 2 emissions, since it compiles all available fluxes, including volcanic lakes, as well as diffuse emissions. We therefore consider our assumptions to be in line with other evidence, although the uncertainties here are high.

Weathering
The role of continental weathering in the glacial-interglacial carbon cycle has been a subject of interest for a long time (see, e.g., Gibbs & Kump, 1994;Munhoven, 2002;Munhoven & François, 1996). Recently, these earlier studies have been revisited by Börker et al. (2020), who also analyzed the role of a loess cover and its changes on the global weathering fluxes. They furthermore used an improved reconstruction for the lithology of the continental shelves exposed at the LGM. Börker et al. (2020) find that neither the global CO 2 consumption nor the global HCO − 3 production rate by continental weathering exhibit significant variations between the LGM, the Mid-Holocene and recent time: all the calculated changes remain within the uncertainty range of each flux. It should nevertheless be noticed that Börker et al. (2020) base their estimates on a continental weathering model, driven by climatic boundary conditions (i.e., runoff) derived from an Earth System Model (ESM). Continental weathering models generally yield comparatively low estimates of global CO 2 consumption and HCO − 3 production rates, even when forced by realistic runoff distribution estimates, compared to estimates based upon river chemistry data: data-based flux estimates are up to 50% greater than those derived from the continental weathering models (see, e.g., the discussion by Munhoven, 2002). In addition, runoff distributions derived from ESMs generally present a low bias. When used directly to force continental weathering models, they therefore contribute to further reduce the global flux estimates. This is also illustrated by Börker et al. (2020) who show that when using the ESM runoff distribution, they obtain a global CO 2 consumption rate of 7.5 Tmol/yr, whereas a river-data based runoff distribution yields 16.8 Tmol/yr; for comparison Gaillardet et al. (1999) find 24.0 Tmol/yr, using a data-based approach. For the concomitant global riverine HCO − 3 (alkalinity) production, Börker et al. (2020) find 11.2 and 23 Tmol/yr, respectively, while Gaillardet et al. (1999) find 36.3 Tmol/yr. For our weathering Scenario W0 (also included in the standard run SE), we adopt here HCO − 3 production rates of 12 Tmol/yr for each of both carbonate and silicate weathering sources (equivalent to a global CO 2 consumption rate of 18 Tmol/yr and a global HCO − 3 production rate of 24 Tmol/yr), which we keep constant in time. For the sake of completeness, and in the light of earlier studies (Gibbs & Kump, 1994;Munhoven, 2002) which suggested that weathering fluxes have not remained constant in time, we nevertheless include four additional scenarios (W1-W4) in our investigation. In each of these, we scale one of the carbonate or silicate weathering related HCO − 3 sources by sea level, so that they vary by ±20% around their mean, with maximum rates either at peak glacial or peak interglacial times (four possible combinations-see Figure 5c). Weathering rates thus vary between 9.6 and 14.4 Tmol/yr, an increase by +50%.
Note that modeling studies focusing on the carbon cycle evolution on time scales of millions of years and more often use CO 2 -dependent weathering flux rate laws, such as W X = W 0 (CO 2 ∕CO 2,0 ) n X , with n Si = 0.2 and n Ca = 0.4 for silicate and carbonate weathering, respectively, where W 0 and CO 2, 0 are the pre-industrial flux rate and CO 2 concentration, respectively (Zeebe, 2012). This would lead for LGM CO 2 concentrations to a reduction in the fluxes by 7% and 13% for silicate and carbonate weathering, respectively. Anomalies of this size would be within the range of our sensitivity scenarios (W1-W4) but have not been chosen in our 10.1029/2020PA004020  (Bereiter et al., 2015). Gray line in (c) is the imposed sea level (Bintanja & van de Wal, 2008). standard setup, since the analyses of Munhoven (2002) and Börker et al. (2020) suggest that processes related to glaciation, such as exposed marine shelf areas, loess, and changes in runoff, are the primary controls on continental weathering during G/IG time scales.

Coral Reefs
Some assumptions on coral reef growth, an important shallow water sink for CaCO 3 , are also necessary in a model that considers the most important carbon fluxes in the late Pleistocene. Indeed, in the so-called coral reef hypothesis (Berger, 1982;Opdyke & Walker, 1992) it had been suggested that the dominant part of the deglacial CO 2 rise might have been caused by the enhanced growth of corals associated with sea level rise. Munhoven and François (1996) assumed coral reef growth to be a function of sea level (setting the extent of the shelf area available) and the rate of sea level change with no growth for sea level rising faster than 15 mm/yr or sinking faster than 5 mm/yr. Furthermore, coral reef growth is a function of the saturation of CaCO 3 in seawater ∼ (Ω CaCO 3 − 1) 1.7 tuned in that way that shallow water accumulation over the last 5 kyr reached 7 ⋅ 10 12 mol CaCO 3 per year, in agreement with Milliman (1993). This approach is used here in Scenario C1. Nonreefal carbonate accumulation and erosion also considered in Munhoven and François (1996) is neglected here, since the resulting fluxes would become too large and in disagreement with other approaches (e.g., Brovkin et al., 2012). The G/IG amplitude during the late Pleistocene in shallow water carbonate accumulation in our Scenario C1 is with ∼8 Tmol/yr (Figure 6c) similar to what has been used in the CLIMBER model (personal communication with authors of Willeit et al., 2019).
Recently, reef carbonate productivity during Quaternary sea level oscillations has been estimated by a numerical model (Husson et al., 2018). When compared to our Scenario C1, it becomes apparent that the different reef carbonate productivity fluxes agree reasonably well along glacial periods but largely differ during terminations (Figure 6c). The accumulation rates calculated by Husson et al. (2018) rise to about 90 Tmol/yr around 420 kyr BP. Coral reefs would thus have accumulated about 1,000 Pmol over a 20 kyr period. Such an extreme large shallow water sink for carbonate is, however, unrealistic since it would consume a large part of the 3,300 Pmol of total alkalinity contained in the ocean (Figure 3c), even if alkalinity supply by continental weathering is taken into account. Taken at face value, coral reef growth based on Husson et al. (2018) cannot be processed with BICYCLE-SE, since the CO 2− 3 concentration in various boxes would around 414 kyr BP decrease to 0, at which point our numerical solution breaks down. We therefore implemented a version of Husson et al. (2018), in which reef productivity is limited to at maximum 20 Tmol/yr (Scenario C2). Note that Husson et al. (2018) did not investigate the consequences of their calculated coral reef growth rates with a carbon cycle model. The authors approximated the effect on atmospheric CO 2 by the rule of thumb that at maximum 0.6 mol of C might be emitted to the atmosphere for each mole of precipitated carbonate referring to Ware et al. (1992).  (Bereiter et al., 2015). Gray line in (c) is the imposed sea level (Bintanja & van de Wal, 2008).

Improvements Due to the Process-Based Sediment Model
In our standard run of model version SE we find the following behavior of the sediment-ocean interaction ( Figure 3): Relevant variables of the water column reproduce the dynamic behavior of the model without sediments (Scenario AOB) remarkable well, although with an offset toward lower numbers. This refers, for example, to the deep ocean CO 2− 3 concentration in both Atlantic and Indo-Pacific, and the total amounts of alkalinity or DIC in the ocean. This agreement between SE and AOB confirms that the simplified expression of the carbonate compensation by a delayed response function in version AOB is a very good approximation of the ocean-sediment interaction. The offset between both scenarios can be understood as follows: In AOB all CaCO 3 that enters the deep ocean is dissolved again, and ocean-sediment fluxes only occur to counteract any offset in CO 2− 3 concentration. In SE all CaCO 3 above (and some below) the saturation horizons is accumulated in the sediments, and thus lost to the ocean, thereby leading to lower total alkalinity and DIC, and consequently also to lower CO 2− 3 . Dynamical changes in the saturation horizon due to any climate change leading to carbonate chemistry variations are in the sediments directly seen in vertical shifts of sediments with the lowest carbonate fractions, for example, the 10% line close to the carbonate compensation depth. Here, small changes in accumulation or dissolution have the strongest impact on this relative measure of carbonate fraction, while the relative changes are smaller at the upper end of the lysocline transition zone, and therefore less obvious in the 50% or 90% lines of carbonate fraction plotted in Figures 3a and 3b. Interestingly, the sediments in the Indo-Pacific hardly record typical glacial/interglacial cycles as has been reconstructed for the equatorial Pacific (Farrell & Prell, 1989). Archer (1991) pointed out that these changes in sediments might be not due to dissolution as suggested by the authors but due to biological production changes. Since marine export production in our model is in the equatorial regions constant over time, no such changes can be recorded in the sediments so far.
We find the following dynamics during glacial terminations, in line with what would be expected from theory (e.g., Sarmiento & Gruber, 2006): After large amounts of CO 2 have been shifted by various processes from the ocean to atmosphere and terrestrial biosphere (similar to a loss of DIC but no loss of alkalinity), oceanic pH is rising, leading to a rise in CO 2− 3 concentration, and therefore to a deepening of CSH by more than 1 km allowing a larger amount of CaCO 3 to accumulate (less to dissolve) above the CSH, leading consequently thereafter to an oceanic loss of CO 2− 3 (equivalent to further losses of DIC and of alkalinity). This dynamic 10.1029/2020PA004020 is also depicted in the carbon contained in the sedimentary mixed layer, which peaks after terminations, synchronous with the drop in deep ocean CO 2− 3 . The changes in the Atlantic basin related to ocean-sediment processes (Figure 3a) are more directly mirroring the G/IG dynamics, probably due to the smaller ocean volume including less amount of DIC and alkalinity, which are more directly affected by already relatively small changes. Some rapid changes in CSH, CO 2− 3 , and the CaCO 3 fraction in the Atlantic appear (e.g., around 10 kyr BP), which are all related to the prescribed instantaneous changes in Atlantic meridional overturning circulation (AMOC) switching from a weak glacial to a strong interglacial mode and vice versa (Figure 2b).
In the deep Atlantic the CO − 3 concentration has its highest values during times of glacial maxima, after which it gradually decreases by 20-30 μmol/kg. In the deep Indo-Pacific the maxima in CO − 3 appears during the deglaciation and at the onset of the interglacials, after which CO − 3 gradually decreases, eventually leading to similar values during glacial maxima and the end of the following interglacial (e.g., LGM vs. pre-industrial). This dynamic leads with 10-20 μmol/kg to a smaller variability of CO − 3 in the deep Indo-Pacific than in the deep Atlantic. Furthermore, during glacial periods CO − 3 in both ocean basins varies on the same order of magnitude than during deglaciations, and again with larger variations in the Atlantic than in the Indo-Pacific.
When we compare the simulated deep ocean CO 2− 3 with reconstructions we have to acknowledge that we can not resolve regional effects, for example, those shown in data from the South China Sea (Luo et al., 2018). However, some general reconstructed CO 2− 3 dynamics from deep ocean sites are worth discussing. So far, deep ocean CO 2− 3 reconstructions are restricted to a few sites from the last 160 kyr (Yu et al., 2013), and one site for the last 250 kyr (Qin et al., 2017) (Figures 3b and 3c). Overall, the pattern of simulated larger changes in CO 2− 3 in the deep Atlantic than in the deep Indo-Pacific basin agrees with reconstructions (Yu et al., 2013(Yu et al., , 2014. The declining values in CO 2− 3 from LGM to pre-industrial times are similarly seen at 1.8 km water depth in the equatorial Atlantic (Yu et al., 2013), while in the very deepest part of the North Atlantic CO 2− 3 seemed to have an increasing trend (Yu et al., 2014). It was shown that to some degree proxy-based CO 2− 3 depends on the foraminifera species chosen for reconstruction (Qin et al., 2017). The last 160 kyr in the Pacific in two proxy records (Qin et al., 2017;Yu et al., 2013) contains similarly small CO 2− 3 variability but in detail different values. While the more shallow core from Qin et al. (2017) showed little change from LGM to preindustrial, there is a similar gradual decline in CO 2− 3 in the deeper core (Yu et al., 2013) during the end of Termination I/beginning of the Holocene than in our simulations. Furthermore, a recent study (Yu et al., 2020) that identified the expansion of carbon-rich Pacific deep waters into the South Atlantic contains in the Antarctic Bottom Water at 5 km water depth a CO 2− 3 signal with peak amplitude of ∼20 μmol/kg and its maximum at the onset of the Holocene, very similar to our simulated deep Indo-Pacific CO 2− 3 . Calculating the CO 2− 3 difference between northern Atlantic and equatorial Pacific as done in Yu et al. (2020) we find (apart from an anomaly around 80-90 kyr BP) a remarkable overlap of our simulations with reconstructions (Figures 3e). In detail, ΔCO 2− 3 contained a G/IG amplitude of ∼40 μmol/kg and is well anti-related to atmospheric CO 2 , suggesting that the most important pattern of changes might indeed been captured with our simple box model. In this comparison with reconstructions one needs to be aware that our box model resolves only changes in pressure with depth, but below 1,000 m water depth DIC and total alkalinity are homogeneous within each ocean basin.  Figures  3c and 3d. Carbon in the mixed sedimentary layer varies on the order of 300 Pg C with peaks during glacial terminations following the described carbonate compensation feedback described above.
The consequence of the implementation of sediments and all the solid Earth fluxes led to in general higher simulated atmospheric CO 2 concentrations (Figure 3e). The regression of modeled versus measured CO 2 is with r 2 = 0.77 in Scenario SE nearly similar than the r 2 = 0.78 obtained in AOB. Thus, the improvement of the model by the implementation of additional processes known to be of relevance in the late Pleistocene leads to hardly any changes in the match of the simulations against the the CO 2 data. However, we have to point out that none of the other processes have been optimized to better fit any data.

CO 2 Release From Volcanism
The effects of the different volcanic emission scenarios on the carbon cycle are as follows (Figure 4): when compared to constant volcanic CO 2 emissions the Scenario V2 with CO 2 release being inversely related to sea level change leads to reductions in atmospheric CO 2 concentration of up to 7 ppm, and to increases of less than 2 ppm. Scenario V2 was motivated by the fact that during times of decreasing obliquity (and therefore of land ice sheet buildup) and on multi-millennial time scales the drop in atmospheric CO 2 lagged the fall in sea level and in global temperature (Hasenclever et al., 2017;Köhler et al., 2018). Accordingly, a process is required that keeps atmospheric CO 2 high while sea level or temperature was falling. This temporal offset between CO 2 and other records during glacial inception has already been established for the last glacial cycle in Barnola et al. (1987). A change of hydrothermal activity in response to sea level change has been suggested to partly explain this offset (Hasenclever et al., 2017). If we now consider this process in its full implication, including its reverse effect during sea level rise, which was ignored by Hasenclever et al. (2017), we find that the idea falls short of explaining the CO 2 -sea level offset during all phases of falling obliquity. Especially during terminations, this process leads to reduced CO 2 emissions when compared to constant volcanism. When glacial inceptions-due to decreasing obliquity-at the end of the interglacials start and sea level begins to drop this process leads to rising CO 2 emissions. However, due to the previous emission reductions during glacial terminations, atmospheric CO 2 has not yet fully recovered to concentrations that would have been obtained with constant volcanic emissions (thus ignoring this sea level dependency).
During certain phases, for example, the MIS 5/4 transition around 80-70 kyr BP, which has been one of the initial cases of the CO 2 -sea level decoupling that motivated the study of Hasenclever et al. (2017), CO 2 in Scenario V2 rises 1 ppm above that in Scenario V0. This is only a fraction of the 18 ppm that needs to be explained here and of the 5-9 ppm that has been proposed by the same model when focusing only on the CO 2 release during sea level fall (Hasenclever et al., 2017). Note that the amount of CO 2 injected by volcanic emissions that stays airborne is in the BICYCLE model at the low end of ranges obtained from a variety of models (Joos et al., 2013), implying that a different model with a more average airborne fraction might produce an anomaly in atmospheric CO 2 of about twice that size, for example, 2 ppm. Furthermore, reconstructions of sea level change, on which the imposed volcanic CO 2 emissions are based on, are also in detail uncertain. For example, in Hasenclever et al. (2017) sea level for the last glacial cycle has been investigated based on a model-based inversion of benthic 18 O (Bintanja et al., 2005, in a similar way also considered here), on corals (Medina-Elizalde, 2013), and on an interpretation of Red Sea data (Grant et al., 2012). For the MIS 5/4 transition these different suggestions have been condensed into four scenarios with rates of sea level change between 60 m in 5 kyr (12 m/kyr) and 80 m in 15 kyr (5 m/kyr). The sea level change imposed here is at the lower end of this range. A combination of a higher airborne fraction and the application of faster rates of sea level change might indeed lead to a higher ΔCO 2 anomaly for the MIS 5/4 transition, of maybe ∼5 ppm. The general problem, however, still exists: reduced volcanic CO 2 emission during rising sea level (terminations) prevents that at glacial inceptions the CO 2 -sea level decoupling is explained by this process.
In other words, atmospheric CO 2 in Scenario V2 most of the time stays below CO 2 in Scenario V0 (or rises at maximum 2 ppm above that in V0), and therefore the idea that this process keeps atmospheric CO 2 concentrations high while sea level falls, is only weakly supported-in certain time windows and with a too small amplitude. This finding, however, does not contradict the overall idea of Hasenclever et al. (2017) that sea level changes influence hydrothermal activity at mid-ocean ridges and hot spot volcanoes. It only implies that this process cannot explain the apparent CO 2 -sea level/temperature offset. Other processes need to be identified to explain this decoupling. A combination of coral reef growth, peat accumulation, and human land use could in three different models (Brovkin et al., 2016) explain CO 2 dynamics including the mentioned CO 2 -sea level/temperature decoupling at the end of the Holocene. These same processes (disregarding anthropogenic land use change) were unable to explain the dynamics at the end of the last interglacials (Brovkin et al., 2016). Therefore, a generic process which explains the observed CO 2 -sea level decoupling during each period of falling obliquity remains elusive.
In Scenario V1, where volcanic CO 2 emission was supposed to be inversely related to continental ice volume changes, we find the following effects with respect to the Scenario V0 that contained a constant volcanic CO 2 emission rate. During terminations and the corresponding reduction in land ice volume the atmospheric 10.1029/2020PA004020 CO 2 concentration rose by more than 20 ppm, while during glacial inceptions (land ice growth) negative CO 2 anomalies are with less than −10 ppm smaller in amplitude. This difference in CO 2 amplitude exists because the land ice sheets grow more slowly than they disintegrate.
Due to the chosen similar response time in both processes the combined effect of land ice change and sea level change on volcanic CO 2 emissions obtained in Scenario V3 is identical to Scenario V1 but with only two thirds of its amplitude: Atmospheric CO 2 (relative to V0) rises by up to ∼13 ppm during terminations. This combined scenario also does not contribute to an explanation of the CO 2 -sea level decoupling seen in the data during falling obliquity. If a different response time would be implemented between sea level change and hydrothermal activity, for example, the 50 kyr assumed in Huybers and Langmuir (2017), the two CO 2 volcanic outgassing fluxes might in the most extreme case add to each other, potentially leading to a CO 2 rise of more than 20 ppm during terminations. However, until further evidences for such response times emerge we refrain from implementing those scenarios in our simulation study.

Weathering
The consequences of the time-dependent weathering rates are symmetric leading to similar but opposite changes in atmospheric CO 2 when weathering rates are either increased or decreased (Figure 5b): The higher (lower) glacial carbonate weathering (W1 and W2) leads to a reduction (increase) of CO 2 by less than 4 ppm with respect to constant weathering rates. For variable silicate weathering (W3 and W4) the amplitudes are with ∼20 ppm by about a factor of 4 higher than for carbonate weathering.

Coral Reefs
In comparison to simulations without coral reef growth the approach applied in C1 with reef growth being a function of sea level change (also used in our standard run SE) leads to an additional rise in atmospheric CO 2 during terminations of 15-20 ppm (Figure 6b). This would imply that coral reef growth is a process contributing about a fifth to the reconstructed G/IG rise in atmospheric CO 2 of 80-100 ppm. Changes in the carbon cycle based on the restricted fluxes following Husson et al. (2018) in our Scenario C2 would lead to G/IG changes of approximately 50 ppm for most terminations and to a maximum of more than 100 ppm around 420 kyr BP. Consequently, all other changes in the carbon cycle would have to sum up to changes in CO 2 of less than 50 ppm during terminations, which seems highly unlikely. We therefore conclude that the suggested fluxes of Husson et al. (2018) cannot be reconciled with other known changes in carbon cycle and seemed to be by an order of magnitude too large.

Model Caveats
Our philosophy here was to combine modules of the carbon cycle which are of a similar level of reduced complexity. Thus, we were not aiming for the most ambitious sediment module that might also consider processes such as the influence of organic carbon or of gradients in the solutes of the carbonate system in the sedimentary mixed layer on CaCO 3 dissolution as implemented in other models (e.g., Brovkin et al., 2007;Heinze et al., 1999;Kurahashi-Nakamura et al., 2020;Ridgwell & Hargreaves, 2007;Tschumi et al., 2011). Our box model approach would then include a lot more details in those sedimentary than in others processes. A similar ocean-sediment exchange scheme than ours was adopted in LOSCAR (Zeebe, 2012).
Furthermore, the constant concentration of all tracers below 1,000 m water depth restricts the ability of our model to capture only first order effects of the sediment-ocean interaction. BICYCLE-SE, however, is more specific and process-based than the so far implemented response function of carbonate compensation to deep ocean CO 2− 3 changes which hindered the model application to questions related to solid Earth processes. In spite of the simplified approach, our simulated changes should roughly be in line with reconstructions. When comparing deep ocean components of the simulated carbon cycle with data one should keep these model limitations in mind and focus on timing and amplitude of temporal changes but not on absolute numbers of the variables. Under these considerations and within the reduced complexity of the approach we find the reasonable agreement of simulate and reconstructed CO 2− 3 as discussed in section 3 quite remarkable.
Apart from these sediment specific details atmospheric CO 2 calculated with box models is known to have a high sensitivity to variations in the high latitudes, while CO 2 derived from more complex Earth system models is more sensitive to changes in the low-latitude oceans. The comparison of BICYCLE with other models has shown that the model sensitivity of the ocean module in BICYCLE, at least of its physical system, is Figure 7. Contribution of solid Earth processes in standard run SE to changes in the global carbon cycle. The relevant process and the carbon species that is changed is mentioned in the legend. All fluxes apart from the volcanic outgassing lead to changes in both carbon in AOB and ocean alkalinity. Net ocean-to-sediment flux (yellow) and sediment accumulation (black) are plotted with a temporal resolution of 2 years (all other records with 100 years resolution) to visualize individual sediment dissolution events that lead to a spiky net ocean-to-sediment CaCO 3 flux (difference between yellow points and black line highlighted by gray area). Volcanic outgassing and coral reef growth are plotted as 1-kyr running means.
in the range of other multibox models and rather close to that of OGCMs. The so-called Harvardton Bear Equilibration Index calculated for a specific comparison of such details (Broecker et al., 1999) is for BICY-CLE in the range typical for OGCMs and higher than that of other multi-box models such as CYCLOPS or PANDORA (see details in Köhler et al., 2005).
Box models have clear limitations when it comes to the response of changes in specific water masses on the carbon cycle. Although, the AOB version of BICYCLE has been used in an effort to understand how temporal changes in northern versus southern deep water formation (e.g., Antarctic Bottom Water vs. NADW) might influence the carbon cycle (Lourantou et al., 2010), such exercises are better performed with OGCMs or full ESMs.

General Discussion
Due to the feedback character of the sedimentary response to any changes in deep ocean carbonate ion concentration and the necessary, on long-term average, balance of volcanic CO 2 degassing and silicate weathering a quantified contribution of these solid Earth processes to changes in atmospheric CO 2 is difficult to calculate. A more solid quantification of the contribution of the different processes to changes in the carbon cycle can be done by the fluxes of carbon, either as CO 2 , HCO − 3 , CO 2− 3 , or CaCO 3 , that are entering of leaving the AOB subsystem of the carbon cycle, as shown in Figures 4c, 5c, and 6c, combined together in Figure 7. The contributions of all four processes-volcanism, weathering, coral reefs growth and ocean-sediment exchange fluxes-are with less than 20 TmolC/y each all of the same order of magnitude. Please remember that carbonate accumulation by coral reefs and volcanic activity are functions of sea level (or sea level change), while weathering is assumed to stay constant, and ocean-sediment interaction is a response or feedback to other changes in the carbon cycle. Here, due to our coarse model resolution rather high annual carbonate fluxes are injected in the ocean for single years following individual sedimentary dissolution events. However, the overall model smooths them out over time and space. They are thus not seen in this spiky form in any other variable.
A compilation of the most important fluxes for preindustrial/present day is given in Table 2, including reconstructed numbers and corresponding references. All-in-all, the applied fluxes fall within the range known from reconstructions. The carbon pool size of the sedimentary mixed layer is with 700 Pg C smaller than the amount CaCO 3 (1,600 Pg C) dissolvable by anthropogenic emissions (Archer et al., 1998). However, in an application of BICYCLE-SE to present and future anthropogenic emissions more than 2,000 Pg C has been released from the dissolution of sedimentary CaCO 3 , pointing to the non-negligible role of the historical sedimentary layers (Köhler, 2020). We furthermore compare our solid Earth fluxes briefly with alternative modeling approaches, namely the most recent carbon cycle changes obtained with the CLIMBER model (Willeit et al., 2019). Here, data from a CLIMBER scenario without trend in volcanic CO 2 outgassing and the best fitting regolith removal scenario are considered. CLIMBER also considers volcanic CO 2 outgassing as a function of land ice volume but time delayed by 5-kyr to sea level change, similar to our Scenario V1 but with smaller magnitude varying between 4 and 8 Tmol/yr. Coral reefs as shallow water sinks vary in CLIMBER between 4 and 14 Tmol/yr, in similar G/IG amplitude as in BICYCLE-SE but offset by 5 Tmol/yr toward higher values. The combined silicate and carbonate weathering flux varies between 33 and 50 Tmol/yr on G/IG timescale in CLIMBER, thus is up to twice as large as the constant 24 Tmol/yr assumed in our W0 scenario. Deep ocean CaCO 3 burial varies between 3 and 21 Tmol/yr in CLIMBER, thus during some time smaller than here but generally in the same range as our simulated output. Altogether, considering the differences in model complexity the solid Earth fluxes of carbon are roughly similar in both approaches.
When compared to other, non-solid Earth processes, it is useful to carry out a factor separation analysis as done before (e.g., Köhler et al., 2010), in which all-but-one processes are active with respect to all processes are active. By that approach the individual contributions from ocean temperature, sea level, sea ice, North Atlantic Deep Water (NADW) formation strength, Southern Ocean deep vertical mixing, iron fertilization of the Southern Ocean marine biology, and terrestrial biosphere have been distinguished. When applied with the extended forcing data set for the last 800 kyr with and without solid Earth processes we obtained more or less identical results as before for the individual contributions ( Figure 8). The only process which generates a response that differs largely from previous applications is sea level change in the SE setup (cyan line in Figure 8d), since now, not only the effect of salinity and ocean volume change are included here but also the sea level-related coral-reef growth and volcanic activity.
Our approach of transient carbon cycle box model simulations is valuable for a rough quantification of the contribution of individual processes to global carbon cycle changes. The most useful information of our results lies in the simulated temporal changes as response to the applied scenarios, not in the absolute numbers of our results. The factor separation approach summarized in Figure 8 has furthermore the difficulty that non-linearities are not well presented. For example, when alternatively one-at-a-time instead of all-but-one processes are investigated the contribution of individual processes to changes in atmospheric CO 2 differs (depending on the process) by a few ppm to more than 10 ppm with ocean temperature being the most sensitive with respect to the chosen factor separation (Köhler & Fischer, 2006). Furthermore, when compared for well defined carbon cycle perturbation experiments (Joos et al., 2013) BICYCLE was found to has an airborne fraction on the lower end of the range span by other, more complex models (Hasenclever et al., 2017), suggesting that especially CO 2 injected into the atmosphere (e.g., during volcanic outgassing) is in BICYCLE taken up by the ocean faster than in other models. In that sense the stated solid Earth contributions to atmospheric CO 2 changes should be understood as an conceptual analysis and the quantified changes as model-specific response to the applied scenarios. Our new standard run SE would lead to a G/IG rise in marine 13 C of 0.53‰ across Termination I-between 23-19 kyr BP (LGM) and 6-0 kyr BP (late Holocene). Peterson et al. (2014) analyzed the difference between the same time windows from a compilation of more than 400 sediment cores and obtained a 13 C rise of 0.38 ± 0.08‰ (2 ) for the ocean between 0.5 and 5.0 km water depth. Making rough assumptions on 13 C in the missing water masses at the surface and in the deep ocean Peterson et al. (2014) finally obtained a whole-ocean change of 0.34 ± 0.19‰. Thus, our results would just meet the upper end of the estimated range. In Scenario SE this change in 13 C would correspond to an increase of the terrestrial carbon storage by 680 Pg C, which is 200 Pg C smaller than the mean change (but within the 1 uncertainty range) suggested by Jeltsch-Thömmes et al. (2019), which has been constrained by simulations using the Bern3D Earth system model and its emulation (including sedimentary response and weathering) to meet certain targets given by various data sets (CO 2 , 13 C in the atmosphere and the ocean, CO 2− 3 ). Parameter optimization within BICYCLE (e.g., of , the parameter responsible for CO 2 fertilization) could be applied in our approach to obtain similarly high terrestrial carbon change, but this would bring marine 13 C further away from the data of Peterson et al. (2014). We therefore refrain from calling upon such optimization approaches.
Recently, the oxidation of pyrite on the continental shelf has been suggested as another process that contributes to a sea level related CO 2 release during the Quaternary (Kölling et al., 2019). For the pyrite oxidation carbonates several tens of meters deep within the sediment are assumed to be used. Such deep sedimentary carbonates are not included in BICYCLE-SE and therefore this process cannot yet be implemented and tested in our model.

Conclusions
In this study, we have implemented processes from the literature to describe solid Earth fluxes within the carbon cycle box model BICYCLE-SE over the past 800 kyr. The model has been extended by a process-based model of early diagenesis in the marine sediments and without further tuning of already implemented processes (marine carbon pumps and terrestrial biosphere) we here investigated the consequences of volcanic CO 2 release, silicate and carbonate weathering, and coral reef growth. We found that previously suggested fluxes of coral reef growth (Husson et al., 2018) are an order of magnitude too large to be in agreement Paleoceanography and Paleoclimatogy 10.1029/2020PA004020 with commonly acknowledged carbon cycle changes. Furthermore, sea level change-related fluxes of volcanic CO 2 release from mid-ocean ridges and hot spot volcanoes can only seldomly and partly explain the multi-millennial decoupling of CO 2 and sea level/temperature as previously suggested (Hasenclever et al., 2017;Köhler et al., 2018). Silicate and carbonate weathering rates are assumed to stay constant over time as recently suggested (Börker et al., 2020) for LGM and modern times; however, this constancy seems in the light of all other changing processes in the carbon cycle rather unlikely. Testing changes in weathering rates by ±20% lead to not yet considered anomalies in atmospheric CO 2 of up to 20 ppm. Our process-based sediment module leads to deep ocean CO 2− 3 dynamics that agree with reconstructions and to quantitatively similar results than the previously applied time-delayed carbonate compensation response function. The latter gives confidence in the results obtained in previous applications of the BICYCLE model without solid Earth processes.

Appendix A: Main Equation on CaCO 3 Dissolution Used in the Sediment Module
The general expression used in the sediment module to quantify the dissolution rate of calcite in the surface seafloor sediment for water depths below the CSH writes where • i = 1, … , 5 denotes the seafloor depth profiles considered (one underneath each surface reservoir); • = 0, … , 79 denotes the sea floor depth intervals [Z j , Z j + 1 ], where Z +1 − Z = 100 m at each time (see Figure 1c). • (R CaCO 3 out ) i, : the flux of carbonate ions in mol/yr out of the mixed layer of a sediment stack (i, j) into the deepest water box lying above; • k sed = 2 ⋅ 10 4 mol m 2 yr ( mol m 3 ) − (n+1) 2 is the dissolution constant; • n = 4.5: dissolution rate order; • a(Z i, j ): is the cumulative seafloor area shallower than Z i, j , in the profile i in m 2 ; • Δ top(Z i, ) = [CO 2− 3 ] sat − [CO 2− 3 ] at water depth Z i, j , where [CO 2− 3 ] sat is the CO 2− 3 concentration at saturation with respect to calcite at depth Z i, j • f i, j : the carbonate fraction in the mixed layer of the sediment stack (i, j); For the depth interval where the CSH is located, that is, for Z CSH such that Δ top (Z CSH ) = 0, and j such that Z j ≤ Z CSH < Z j + 1 , the Equation A1 is adapted to read (R CaCO 3 out The details about how Equation A1 was derived from the general early diagenesis equation (including the fundamental hypothesis, the adopted simplifications and approximation, the calibration, and a comprehensive discussion and comparison with similar schemes) can be found in Munhoven (1997, pp. 123-142).