Modeling the evolution of pulse-like perturbations in atmospheric carbon and carbon isotopes: the role of weathering–sedimentation imbalances

. Measurements of carbon isotope variations in climate archives and isotope-enabled climate modeling advance the understanding of the carbon cycle. Perturbations in atmospheric CO 2 and in its isotopic ratios ( δ 13 C, (cid:49) 14 C) are removed by different processes acting on different timescales. We investigate these differences on timescales of up to 100 000 years in pulse-release experiments with the Bern3D-LPX Earth system model of intermediate complexity and by analytical solutions from a box model. On timescales from years to many centuries, the atmospheric perturbations in CO 2 and δ 13 CO 2 are reduced by air–sea gas exchange, physical transport from the surface to the deep ocean, and by the land biosphere. Isotopic perturbations are initially removed much faster from the atmosphere than perturbations in CO 2 as explained by aquatic carbonate chemistry. On multimillennial timescales, the CO 2 perturbation is removed by carbonate compensation and silicate rock weathering. In contrast, the δ 13 C perturbation is removed by the relentless ﬂux of organic and calcium carbonate particles buried in sediments. The associated removal rate is signiﬁcantly modiﬁed by spatial δ 13 C gradients within the ocean, inﬂuencing the isotopic perturbation of the burial ﬂux. Space-time variations in ocean δ 13 C perturbations are captured by principal components and empirical orthogonal functions. Analytical impulse response functions for atmospheric CO 2 and δ 13 CO 2 are provided. Results suggest that changes in terrestrial carbon storage were not the sole cause for the abrupt, centennial-scale CO 2 and δ 13 CO 2 variations recorded in ice during Heinrich stadials HS1 and HS4, though model and data uncertainties prevent a ﬁrm conclusion. The δ 13 C offset between the Penultimate Glacial Maximum and Last Glacial Maximum reconstructed for the ocean and atmosphere is most likely caused by imbalances between weathering, volcanism, and burial ﬂuxes. Our study highlights the importance of isotopic ﬂuxes connected to weathering–sedimentation imbalances, which so far have been often neglected on glacial– interglacial timescales.


Introduction
Research efforts to understand the past and future evolution of atmospheric CO 2 and of its isotopic signature have been ongoing for many decades. Yet, many open questions remain, including, for example, the offsets in the stable isotope signature of CO 2 (δ 13 C atm ) and of dissolved inorganic carbon (δ 13 C DIC ) between the Penultimate Glacial Maximum (PGM) and Last Glacial Maximum (LGM), as well as between the last (Eem) and current (Holocene)interglacials (Schneider et al., 2013;Eggleston et al., 2016;Hoogakker et al., 2006;Oliver et al., 2010). Could this PGM-LGM δ 13 C offset result from a change in terrestrial carbon storage or rather from changes in carbon and carbon isotope fluxes connected to the sedimentary weathering-burial cycle? Further questions remain regarding CO 2 and δ 13 C atm variations on centennial and millennial timescales and a possible role of changes in the land biosphere carbon inventory (e.g., Köhler et al., 2010;Skinner et al., 2010Skinner et al., , 2017Schmitt et al., 2012;Bauska et al., 2018). Reconstructed variations in CO 2 and δ 13 C atm during Heinrich stadials 4 and 1 (HS4 and HS1), two Northern Hemisphere cold periods during the last glacial, have been attributed to changes in terrestrial carbon storage Published by Copernicus Publications on behalf of the European Geosciences Union.

424
A. Jeltsch-Thömmes and F. Joos: Evolution of atmospheric carbon and δ 13 C perturbations (e.g., Bauska et al., 2016Bauska et al., , 2018; however, other scenarios might be possible. The perturbations in CO 2 and δ 13 CO 2 due to a carbon input (or removal) from the land biosphere or from fossil fuel burning are eventually removed by a few principal processes. In the case of CO 2 , over the first decades and centuries, this removal is because of ocean and land uptake; on millennial timescales, it is CaCO 3 compensation; finally, on timescales of hundreds of thousands of years, the remaining atmospheric CO 2 perturbation is removed by enhanced silicate rock weathering (e.g., Archer et al., 1998;Joos et al., 2001Joos et al., , 2013Colbourn et al., 2015;Lord et al., 2016). In the case of δ 13 C, the perturbation is mixed rather quickly throughout the atmosphere-land-ocean reservoir and is removed on multimillennial timescales by changes in the isotopic signature of burial fluxes to the lithosphere.
The simulation of climate and CO 2 over thousands of years with dynamic three-dimensional models remains challenging. Over the last decades, computational power has continuously increased and first results from simulations with Earth system models of intermediate complexity (EMICs) covering glacial-interglacial cycles have been published (e.g., Brovkin et al., 2012;Menviel et al., 2012;Ganopolski and Brovkin, 2017). However, to date these long timescales prevented large ensemble simulations with spatially resolved models as well as glacial-interglacial simulations using state-of-the-art Earth system models. A possible tool to overcome these computational obstacles is the use of reducedform models invoking impulse response functions (IRFs) and principal component and empirical orthogonal function (PC-EOF) analysis. IRFs and PC-EOF analysis are instrumental to characterize the fundamental time and spatial scales of model response to a perturbation.
The response of a linear system to a perturbation such as carbon emissions to the atmosphere can be fully characterized by its IRF. The atmosphere-ocean CO 2 cycle can be approximated as a nearly linear system, as long as concentrations do not vary much, i.e., less than a doubling of preindustrial CO 2 concentrations . The perturbation in the atmosphere or ocean is then given by the convolution integral of the emission history and the corresponding IRF (e.g., . IRFs can easily be described by a sum of exponentials, yielding a cost-efficient substitute model for describing the temporal evolution of an atmospheric CO 2 perturbation. IRFs have been widely used to calculate remaining atmospheric fractions of anthropogenic CO 2 emissions on different timescales and to characterize model behavior and their characteristic response timescales (e.g., Siegenthaler and Oeschger, 1978;Maier-Reimer and Hasselmann, 1987;Sarmiento et al., 1992;Siegenthaler and Joos, 1992;Enting et al., 1994;Archer et al., 1997Archer et al., , 1998Archer et al., , 2009Ridgwell and Hargreaves, 2007;Archer and Brovkin, 2008;Joos et al., 2013;Colbourn et al., 2015;Lord et al., 2016). However, peer-reviewed studies that quantify the IRF for carbon iso-topes of CO 2 are scarce (Schimel et al., 1996; and, to our knowledge, not available for multimillennial timescales. Principal component analysis is widely used to capture the information of multidimensional data. The atmospheric perturbation in CO 2 and its isotopes in response to carbon input external to the ocean-atmosphere system is propagated to the ocean. The resulting spatiotemporal (4-D) evolution of marine tracers may be conveniently approximated using PC-EOF analysis. The perturbation is described as a superposition of spatial empirical orthogonal functions (EOF(x)) multiplied with the corresponding time-dependent principal component (PC(t)). Combining IRF and PC-EOF then allows one to compute the perturbation in space and time with minimal computational resources Joos et al., 2001) The main goal of this study is to investigate the timescales and underlying processes for the removal of an atmospheric perturbation in CO 2 and its δ 13 C signature in response to an external carbon input. To this end, we performed pulserelease (uptake) experiments over 100 000 years with the Bern3D-LPX EMIC. We determine characteristic impulse response (or Green's) functions (IRF) and describe marine spatiotemporal variations using PC-EOF analysis. Responses of the atmospheric radiocarbon signature ( 14 C) are briefly addressed. Implications of our results for the PGM-LGM δ 13 C difference and for past centennial CO 2 variations are discussed. Simple expressions that illuminate the fundamental differences in the removal timescales for an atmospheric perturbation in CO 2 versus the associated perturbation in δ 13 CO 2 are developed in the Appendix. The pulse experiments represent (i) the sudden carbon uptake from or release to the atmosphere by the land biosphere, (ii) the release of old organic carbon, e.g., from ancient peat or permafrost material, (iii) the burning of all conventional fossil fuel resources, and (iv) a variation in 14 C production by cosmic radiation or atomic bomb tests.

The Bern3D-LPX model
For this study, the Bern3D v2.0s Earth system model of intermediate complexity (EMIC) is used coupled to the dynamic vegetation model LPX-Bern v1.4. The Bern3D model consists of a three-dimensional geostrophic ocean (Edwards et al., 1998;Müller et al., 2006) with an isopycnal diffusion scheme and Gent-McWilliams parameterization for eddy-induced transport (Griffies, 1998) and a single-layer energy-moisture balance atmosphere to which a thermodynamic sea-ice component is coupled (Ritz et al., 2011). We use the ocean model coupled to a 10 layer sediment module (Tschumi et al., 2011;Heinze et al., 1999) with updated sediment parameters (Jeltsch-Thömmes et al., 2019) and to a global representation of rock weathering (Colbourn et al., A. Jeltsch-Thömmes and F. Joos: Evolution of atmospheric carbon and δ 13 C perturbations 425 2013). The resolution of the Bern3D model components is 41 × 40 horizontal grid cells with 32 logarithmically spaced depth layers in the ocean. Wind stress at the ocean surface is prescribed from the NCEP/NCAR monthly wind stress climatology (Kalnay et al., 1996), and air-sea gas exchange and carbonate chemistry follow OCMIP-2 protocols (Najjar and Orr et al., 1999Orr et al., , 2017Wanninkhof, 2014;Orr and Epitalon, 2015). The global-mean air-sea gas transfer rate is reduced by 19 % to match observational radiocarbon estimates (Müller et al., 2006), and gas transfer velocity scales linearly with wind speed (Krakauer et al., 2006).
The sediment module dynamically calculates the transport, redissolution or remineralization, and bioturbation of solid material, the porewater chemistry, and diffusion in the top 10 cm of the sediment. The governing equations of the sediment model are given in Appendix A of Tschumi et al. (2011). A comparison of simulated versus observation-based sediment composition for the setup used in this study is given in Appendix B of Jeltsch-Thömmes et al. (2019). A table with all parameters applied here and by Jeltsch-Thömmes et al. (2019) is given in the Appendix (Table B1). Four solid tracers (CaCO 3 , opal, particulate organic carbon (POC), clay) and seven tracers in the porewater (DIC, DI 13 C, DI 14 C, alkalinity, phosphate, oxygen, and silicic acid) are modeled. The dissolution of CaCO 3 and POC depends on the respective weight fraction of CaCO 3 and POC in the solid phase of the sediment and the porewater CO 2− 3 and O 2 concentrations, respectively. 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;Jeltsch-Thömmes et al., 2019). The model assumes conservation of volume, i.e., the entire column of the sediments is pushed downward if deposition exceeds redissolution into porewaters. 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). The burial flux at preindustrial steady state is 0.22 GtC yr −1 in the form of CaCO 3 , 0.24 GtC yr −1 in the form of POC, and 6.96 Tmol Si yr −1 in the form of opal, all within the observational range. The global inventories in the interactive sediment layers amount to 916 GtC for CaCO 3 , to 510 GtC for POC, and to 21 460 Tmol Si for opal. Carbonate chemistry within sediment porewaters is calculated as in the ocean by using the MOCSY routine of Orr and Epitalon (2015).
Input fluxes by "weathering" of P, ALK, DIC, DI 13 C, and Si are added uniformly to the coastal surface ocean. At the beginning of transient simulations, the global input fluxes of these compounds are set equal to the burial fluxes diagnosed at the end of the model spin-up. These input fluxes are jointly denoted as "weathering flux". These fluxes are further attributed to weathering of organic material, of CaCO 3 , of CaSiO 3 , and to volcanic CO 2 outgassing. The flux of phosphorus (P) is assigned to weathering of organic material, and related carbon (C) and ALK fluxes are computed by multiplication of the P flux with the Redfield ratios where P : C : Alk is equal to 1 : 117 : −17 for organic material. Similarly, the silicon (Si) flux is assigned to CaSiO 3 weathering, and the related ALK flux is computed using Si : ALK is equal to 1 : 2 based on the simplified equation for CaSiO 3 weathering: (Colbourn et al., 2013). The remaining ALK flux is attributed to CaCO 3 weathering with the stoichiometric ratio C : ALK equal to 1 : 2 following from The volcanic flux is the remaining flux needed to balance the C input flux. The diagnosed fluxes at the end of the spin-up are 0.24 GtC yr −1 for weathering of organic material, 0.14 GtC yr −1 for CaCO 3 weathering, 0.08 GtC yr −1 for volcanic CO 2 outgassing, and 6.96 Tmol Si yr −1 for CaSiO 3 weathering. The isotopic signature of the weathering carbon corresponds to the respective signature of the burial fluxes and amounts to δ 13 C = −9.2 ‰, intermediate between isotopically light organic carbon (δ 13 C = −20.5 ‰) and heavier CaCO 3 (δ 13 C = 2.9 ‰). During experiments, weathering fluxes of CaCO 3 and CaSiO 3 , and accordingly DIC, DI 13 C, and alkalinity, are allowed to vary as a function of surface air temperature, runoff, and net primary productivity (NPP). The parameterization of this weathering feedback follows equations described in detail in Colbourn et al. (2013). The runoff dependence is parameterized as a function of temperature, and for the productivity feedback the NPP from the coupled LPX-Bern model is used. We apply weathering feedbacks in the global average 0-D version (Colbourn et al., 2013(Colbourn et al., , 2015; the equations are given in Appendix B. The Bern3D model is coupled to the Land surface Processes and eXchanges (LPX) model v1.4 (Lienert and Joos, 2018), which is based on the Lund-Potsdam-Jena (LPJ) model (Sitch et al., 2003). In LPX, coupled nitrogen, water, and carbon cycles are simulated and vegetation composition is determined dynamically and represented with 20 (15 on natural and 5 on anthropogenically used land) plant functional types (PFTs). The PFTs compete within 426 A. Jeltsch-Thömmes and F. Joos: Evolution of atmospheric carbon and δ 13 C perturbations their bioclimatic limits for resources. CO 2 assimilation of plants is implemented following Farquhar et al. (1980) and Haxeltine and Prentice (1996a, b), and isotopic discrimination during photosynthesis is calculated according to Lloyd and Farquhar (1994). The carbon cycles of Bern3D and LPX-Bern are coupled through carbon and isotope exchange fluxes between the land and atmosphere and between the ocean and atmosphere. Climate change information from the Bern3D model is passed on to LPX-Bern via a patternscaling approach (Stocker et al., 2013). Spatial anomaly patterns in monthly temperature and precipitation, derived from a 21st century simulation with the Community Climate System Model (CCSM4), are scaled by the anomaly in global monthly mean surface air temperature as computed interactively by the Bern3D. These anomaly fields are added to the monthly baseline (1901 to 1931 CE) climatology of the Climate Research Unit (CRU) (Harris et al., 2014). A more detailed description of the LPX-Bern model can be found in Keller et al. (2017) and Lienert and Joos (2018).

Experimental protocol
The Bern3D model is spun-up over 60 000 years (1000 years is represented by "kyr") to a preindustrial steady state with 1850 CE boundary conditions. Atmospheric CO 2 is prescribed to 284.7 ppm, δ 13 C to −6.305 ‰, and 14 C to 0 ‰. For LPX, the land area under anthropogenic use is fixed to its 1850 CE extent and kept constant in all simulations. The LPX-Bern model had an uncoupled spin-up under identical 1850 CE boundary conditions over 2500 years. Bern3D and LPX-Bern are then coupled and equilibrated for 500 years, again under fixed 1850 CE boundary conditions. After the spin-up, atmospheric CO 2 and its isotopic ratio are calculated interactively, with enabled carbon-climate feedbacks. Pulse experiments are started from the 1850 CE equilibrium state at nominal year −100. After 100 years, i.e., during nominal year 0 of the simulations, an external flux of carbon with characteristic isotopic signatures is removed from or added to the model atmosphere at a constant rate over the year. Runs are continued for 100 kyr to simulate the redistribution of the added carbon and isotopes in the Earth system. Table 1 summarizes run names and key characteristics of all simulations. Pulse sizes are varied between −250 GtC (removal) to +500 GtC (addition) with δ 13 C of −24 ‰ and 14 C of 0 ‰. This corresponds to a hypothetical, sudden uptake or release of carbon from the land biosphere with a typical C 3 -like δ 13 C signature. For simplicity, and to ease interpretation, we set the 14 C of these land biosphere carbon pulses to the atmospheric signature; thus, we assume that the released material has been recently assimilated and neglect the small depletion in 14 C of "young" plant, litter, and soil material. In sensitivity experiments, 14 C is set to −500 ‰ and −1000 ‰ to mimic the release of old, dead, organic carbon, e.g., from buried peat or permafrost soils instead of young plant-derived material. In another simulation, the iso- topic signatures are set to the signature for the modern mix of fossil fuels. Finally, we also performed a run where only 14 C, but no 12 C or 13 C, is added to represent a change in radiocarbon production in the atmosphere. All these simulations were performed with fully interactive ocean sediments and enabled CaCO 3 and CaSiO 3 weathering feedbacks. The influence of the weathering feedback and of ocean-sediment interactions is quantified using two additional simulations. First, the weathering input flux is kept time invariant ("no weathering feedback"), Second, the sediment module is, in addition, not included in a so-called "closed" atmosphereocean-land model setup.
We also assess the influence of glacial climate boundary conditions and of a different land biosphere model in two additional simulations: LGM 500 and 4box 500 . In 4box 500 the Bern3D setup and spin-up is as in WEA 500 , except that an atmospheric CO 2 concentration of 278 ppm is used instead of 284.7 ppm, and the pulse is released in the first time step of simulation year 100, rather than being distributed over year 100. The model configuration includes reactive sediments and weathering feedbacks. A four-box terrestrial biosphere model (Siegenthaler and Oeschger, 1987), that allows for CO 2 fertilization, instead of LPX-Bern is coupled to the Bern3D in these runs. In LGM 500 the model is forced by glacial boundary conditions instead of preindustrial conditions as in 4box 500 and in WEA 500 . CO 2 is set to 180 ppm and Northern Hemisphere ice sheets are set to Last Glacial Maximum coverage (Peltier, 1994). It has to be noted, however, that forcing the model with 180 ppm during spin-up leads to less carbon stored in the ocean under LGM compared to preindustrial (PI) conditions, and results from LGM 500 should be treated with some caution. Differences in results between WEA 500 and 4box 500 are due to differences in the two land biosphere models. Differences between 4box 500 and LGM 500 are exclusively due to differences in climatic boundary conditions.

Data reduction
The remaining fraction of a pulse-like perturbation or IRF at time t after the pulse release is defined for atmospheric CO 2 as indicates a perturbation, here evaluated as the difference between a simulation with pulse release and a corresponding control simulation. P indicates the pulse size in gigatonnes of carbon (GtC), which yields the initial atmospheric CO 2 perturbation in parts per million (ppm) when divided by 2.12 GtC ppm −1 . The superscript "ctrl" refers to the control simulation without pulse, "ini" to the (maximum) initial perturbation assuming an instantaneous carbon release, and subscript "a" to the atmosphere. The IRF for δ 13 C a is IRF(δ 13 C a )(t) = δ 13 C a (t) − δ 13 C ctrl a (t) δ 13 C ini a .
(2) δ 13 C ini a denotes the initial isotopic perturbation in the atmosphere and equals: δ 13 C P indicates the isotopic signature of the pulse, the subscript 0 the value prior to the pulse, and N a the atmospheric carbon inventory. Analogous equations to Eqs.
(2) and (3) hold for 14 C. An analytical expression for the IRF is calculated from experiment WEA 500 . The IRF is fitted by the sum of five exponential terms to using a least-squares optimization routine in Python. The coefficients a i are fractions of the perturbation, each associated with a timescale τ i , and a 0 is the fraction of the perturbation remaining constant in the atmosphere (τ 0 = ∞). The sum over all six coefficients, a i , equals 1. Modeled values are well reproduced with n = 5 exponents. Differences between the original model output and the fit are insignificant for practical applications. Finally, we note that the IRF or remaining atmospheric fraction is in our experimental setup always smaller than one: the carbon is added at a constant rate during year 0 in the pulse simulations, and part of this carbon is already taken up by the ocean during this initial year. Similarly, the IRF for δ 13 C a and 14 C a is fitted with n = 4 exponents. Principal component empirical orthogonal function (PC-EOF) analysis is used to fit the 13 C isotopic perturbations of dissolved inorganic carbon (DIC) in the ocean ( δ 13 C DIC ). To this end, the Python package eofs (Dawson, 2016) is used to calculate principal component time series (PC(t)) and the spatial empirical orthogonal functions (EOF(x)) from the volume-weighted 4-D field of δ 13 C DIC (t, x). PCs and EOFs are then used to build a cost-efficient substitute model of the spatiotemporal evolution of δ 13 C DIC in response to an atmospheric pulse perturbation. The output frequency for marine 3-D tracer fields after the perturbation is every 10 years during the first 1 kyr, every 200 years until 10 kyr, and every 1 kyr thereafter. In the PC-EOF calculation this gives more weight to the first part of the simulation, where changes in the perturbation are larger.

Results
3.1 Earth system response to a pulse release of carbon Impulse response functions for atmospheric CO 2 and δ 13 C We first describe the IRF for the atmospheric CO 2 perturbation using experiment WEA 250 (Fig. 1a). The initial spike in response to the carbon release is followed by a substantial decline over the first decades due to carbon uptake by the upper ocean and land biosphere. The IRF is reduced to less than 40 % within the first century and decreases further on the timescale of ocean mixing to less than 16 % by year 2000. The fit by exponentials ( Table 2) yields timescales of 6, 47, and 362 years for this initial decline. They represent the continuum of carbon overturning timescales within the ocean and the land biosphere as well as timescales governing air-sea and air-land carbon exchanges.
CaCO 3 compensation and weathering-burial imbalances remove the atmospheric perturbation on multimillennial timescales such that after 10 kyr less than 8 % and after 100 kyr less than 2 % of the CO 2 perturbation remain airborne (Fig. 1a). The fit by exponentials shows that a fraction of 11.2 % of the initial perturbation is removed with a  Table 1 and Sect. 2.2 for details of the experimental setup of the runs. Data are filtered with a moving average of 1000 years in the interval 1-10 kyr and 5000 years in the interval 10-100 kyr for better visibility. Factorial simulations (WEA 500 , SED 500 , CLO 500 , 4box 500 , and LGM 500 ) for the time interval 1 to 100 kyr are shown in Fig. 2 for better visibility. timescale of 5.5 kyr. We associate these values with the process of CaCO 3 compensation. A remaining fraction of 4 % to 5 % is removed with a timescale of 69 kyr (Table 2), linked to enhanced silicate rock weathering. The jumps visible at around 20 kyr in Fig. 1 arise from a sea-ice-albedo feedback in the control run. They appear most pronounced in the experiments with small perturbations. The change in the control run in CO 2 is less than 1 ppm and less than 0.03 ‰ in δ 13 CO 2 . These jumps are not of further relevance for our discussion. In comparison to others, our results show, up to a few percent, lower remaining fractions for the 5000 GtC pulse (WEA 5000 ) than results presented by Colbourn et al. (2015) and Lord et al. (2016). The remaining airborne frac-tion of the CO 2 perturbation is 5 % and less after 100 kyr in all studies.
Compared to CO 2 , the atmospheric δ 13 C perturbation is initially removed much faster (Fig. 1b). Within the first decades, IRF(δ 13 C a ) is reduced to below 10 %. The further decline is slower, and 1 kyr after the pulse less than 3 % of the initial perturbation remains airborne. The remaining atmospheric perturbation is finally reduced to ≤ 1 % of the initial perturbation at the end of the simulation (Fig. 1b). The fit by exponentials for IRF(δ 13 C a ) ( Table 2) shows that 84 % of the perturbation are removed with a timescale of 6 years, reflecting the initial fast decline seen in Fig. 1b when the perturbation is taken up by the upper ocean and living vegetation. About 9 % and 3 % are removed on timescales of 75 and 436 years, respectively, reflecting the decadal and centennial uptake of the perturbation by the land biosphere and deeper ocean. Approximately 3 % of the perturbations are finally removed with a timescale of ∼ 74 kyr by weathering-sedimentation processes, in agreement with the e-folding timescale of the decline (∼ 70 kyr; see further below). However, it has to be kept in mind that these removal timescales result from a fit, and thus other equally accurate solutions are possible to approximate the continuum of removal timescales.

Processes
Next, we address the processes behind the different temporal evolution of CO 2,a and δ 13 C a . We focus on the ocean as ocean uptake is thought to dominate on long timescales. A perturbation in atmospheric CO 2 or in 13 CO 2 is communicated by air-sea gas exchange to the surface ocean and further transported, mainly by physical processes, to the deep ocean. The fundamental difference between the IRF for CO 2,a and δ 13 C a is linked to the aquatic carbonate chemistry and the associated equilibrium between dissolved CO 2 and bicarbonate and carbonate ions ( son et al., 2007). The removal of an atmospheric CO 2 perturbation is co-controlled by this acid-base carbonate chemistry (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, because 13 CO 2 and 12 CO 2 are affected approximately equally by the acid-base chemical reactions. The chemical buffering diminishes the uptake capacity of the ocean for CO 2 but not for the isotopic ratio. Correspondingly, the IRF for CO 2,a decreases less rapidly than the IRF for δ 13 C a , and IRF(δ 13 C a ) is smaller than IRF(CO 2,a ) over the simulation period.
For an illustrative, semiquantitative analysis, we developed approximate expressions for the IRF for CO 2,a and δ 13 C a for timescales up to 2 kyr as explained in the Appendix. Simplifying assumptions are that equilibrium between the atmosphere and a fraction of the ocean with the initial carbon inventory N o,0 is assumed, while land biosphere carbon uptake and other climate-carbon-cycle feedbacks are neglected. Changes in land biosphere carbon storage can be considerable, in particular during the first century (see Fig. 3a). However, the role of the changes in land carbon stock becomes smaller when more carbon is taken up by the ocean, and thus for simplicity it is not considered in the analytical expression for CO 2 . For CO 2 , the expression reads as The subscripts "a" and "o" indicate the atmosphere and ocean, respectively; 0 represents the time prior to the pulse; and ∞ represents the time when an effective ocean volume with the initial carbon inventory of N o,0 is equilibrated with the atmosphere. ξ describes the influence of the acid-base carbonate chemistry on the relationship between the relative perturbation in the CO 2 partial pressure and in dissolved inorganic carbon (DIC) (Revelle and Suess, 1957). ξ varies with environmental conditions and increases with the size of the CO 2 perturbation. ξ is on the order of 10 for small pulse sizes (here −250 to 500 GtC). The dependency of IRF(CO 2 ) on the size of the carbon pulse is implicitly captured in the variable ξ . The key information of Eq. (5) is that the magnitude of the IRF is determined by the ratio of the (equilibrated) ocean carbon inventory and the buffer factor (N o,0 /ξ ). The smaller this ratio, the larger IRF(CO 2,a ) is.
For δ 13 C a , we also consider the initial carbon inventory of the land biosphere (N b,0 ), as a substantial amount of the isotopic perturbation is contained in the land biosphere on multicentennial-to-millennial timescales (Fig. 3b). The corresponding expression reads as follows (for the derivation, see the Appendix): Equation (6) shows that IRF(δ 13 C a ) is described by the product of two ratios. The first ratio corresponds to the ratio of the carbon inventory in the atmosphere immediately after the pulse release to the total carbon inventory in the ("equilibrated") atmosphere-ocean-land-biosphere system. It tells us that the pulse perturbation is mixed uniformly (see also Appendix Eq. A20) and in proportion to the carbon inventories of the different reservoirs (atmosphere, ocean, land biosphere). The second fraction is about 1.2 and is thus a correction term of order 20 % for small pulse sizes, as typically applied in this study. It arises as the isotopic signature of the ocean and land biosphere is different from that of the atmosphere. δ 13 C mean corresponds to the weighted mean signature of the atmosphere, ocean, and land biosphere, with the fraction of the pulse found in each of these reservoirs as weights (Appendix Eq. A27). Remarkably, the buffer factor ξ does not directly show up in Eq. (6). The reason is, as noted already above, that CO 2 and 13 CO 2 are both about equally affected by the aquatic chemistry. In turn, their ratio, i.e., δ 13 C, is hardly affected by this acid-base buffering. Equations (5) and (6) reveal the fundamental difference between a perturbation in CO 2 (and 13 CO 2 ) and in the isotopic ratio (δ 13 CO 2 ). For small pulse sizes (P → 0), neglecting the influence of the land biosphere and the correction term in Eq. (6) for δ 13 C, Eqs. (6) and (5) are formally equal, except that a buffer factor ξ of 1, instead of ∼ 12, applies for δ 13 C a . Equations (6) and (5) show that the pulse perturbation is distributed between the atmosphere and ocean proportionally to their initial carbon inventories N a,0 and N o,0 Table 2. Fit parameters for IRF( CO 2 (t)) and IRF( δ 13 C a (t)) of WEA 500 , IRF( 14 C a,dead (t)), and IRF( 14 C a,only (t)). Note that in the case of IRF( 14 C a,dead (t)) the actual fitting was done only with data until simulation year 15 000 to avoid complications from the sea-ice feedback in the control run (see main text and Fig. 1c). in the case of δ 13 C a but proportionally to N a,0 and N o,0 /ξ in the case of CO 2 . As a consequence, the perturbation in the ratio δ 13 C a is apparently removed much faster than the perturbation in the concentration CO 2,a , despite the fact that the timescales to transport a perturbation from the surface ocean to the deep ocean by advection, convection, or diffusion are the same. We may illustrate this difference in the removal rate numerically. ξ is about 12 for a small pulse, N a,0 is 600 GtC, and N o,0 is 37 400 GtC for the whole ocean. Further, the top 300 m with a carbon inventory of 2600 GtC is ventilated within approximately a decade through air-sea gas exchange and upper-ocean physical transport processes (circulation, mixing). Then, the evaluation of Eq. (5) means, for a small pulse, that around 73 % and 16 % of the initial CO 2 perturbation are still found in the atmosphere after about 10 years and 2 kyr, respectively. It has to be noted here that we assume ξ to be constant for small perturbations. With increasing pulse sizes, ξ also increases, leaving larger airborne fractions (cf. WEA 500 versus WEA 5000 in Fig. 1a after 1-2 kyr). Regarding δ 13 C, the corresponding fractions (Eq. 6) are much smaller (12 % to 31 % and 1 % to 3 %) for small pulse sizes (−250 to 500 GtC; see Appendix text after Eq. A14 for more details). We assume the perturbation to have mixed with 500 and 2700 GtC on land after 10 yr and 2 kyr, respectively. These values are roughly in agreement with the values shown in Fig. 1a for pulse sizes of up to 500 GtC.
Turning to intermediate timescales (2 kyr ≤ t ≤ 20 kyr), CaCO 3 compensation further reduces the atmospheric CO 2 perturbation and removes the differences in IRF(CO 2,a ) arising from different pulse sizes (Fig. 1a). The process of CaCO 3 compensation is briefly explained as follows. CO 2 is taken up by the ocean and partly reacts to form bicarbonate (HCO − 3 ) and hydrogen (H + ) ions. This makes the water more acidic; the carbonate ion concentration ([CO 2− 3 ]) is reduced and the water becomes corrosive to CaCO 3 . In turn, sedimentary CaCO 3 dissolves to CO 2− 3 and Ca 2+ ions, partly removing the perturbations in [CO 2− 3 ] and [H + ], a so-called seafloor neutralization. Under more acidic conditions, less alkalinity and carbon, in the form of CaCO 3 , is buried in the lithosphere than delivered by weathering input; this excess input leads to an increase in alkalinity and DIC, further removing the perturbations in [CO 2− 3 ] and [H + ] , a socalled terrestrial neutralization, and a new equilibrium between weathering and burial fluxes is established. This chain of processes raises the alkalinity twice as much as the DIC concentration in the ocean, and additional CO 2 is taken up from the atmosphere.
Again, an approximate expression for the IRF(CO 2,a ) is developed. This applies for the time when CaCO 3 compensation for the pulse perturbation is completed (see Appendix): The factor 2 in Eq. (7) arises from the acid-base carbonate chemistry and the factor 1.5 from the change in the ocean carbon inventory due to excess dissolution (or burial) of CaCO 3 , equivalent to ∼ 0.5 × P . The factor 1.5 is model dependent (see Archer et al., 1998) and also slightly different from previous Bern3D model versions. This equation shows that the airborne fraction after CaCO 3 compensation depends on the ratio of the initial atmosphere-to-ocean carbon inventory (N a,0 to N o,0 ). Therefore, IRF CaCO 3 is (approximately) independent of the pulse size. Equation (7) successfully explains why IRF(CO 2,a ) for pulse sizes from −250 to +5000 GtC take on a nearly identical value of 5 % around 20 kyr after the pulse. IRF(δ 13 C a ) stays nearly constant between 1 and 10 kyr (Fig. 1b). This implies that CaCO 3 compensation has little influence on the removal of the δ 13 C perturbation. Reasons are that the change in ocean carbon inventory due to CaCO 3 dissolution (burial) is small compared to the total ocean inventory, and that the isotopic signature of CaCO 3 is similar to that of DIC.
On even longer timescales, the CaSiO 3 weathering feedback further removes the remaining perturbation in CO 2,a with an e-folding timescale on the order of 70 kyr. The perturbation in δ 13 C a is removed by the burial flux of organic LGM 500 . Shown is the time interval 1 to 100 kyr. Data are filtered with a moving average of 1000 years in the interval 1-10 kyr and 5000 years in the interval 10-100 kyr for better visibility. Note that experiment CLO 500 was run for 90 kyr instead of 100 kyr. carbon and CaCO 3 . This burial flux removes the isotopic perturbation to the lithosphere. Surprisingly, the removal rate is substantially faster than expected from the residence time of carbon, and it is highly sensitive to the gradients in the perturbation of δ 13 C ( δ 13 C) within the ocean as further discussed below in the context of the factorial simulations and in the Appendix.

Factorial simulations
We now turn to the factorial simulations to further quantify the effect of CaCO 3 compensation, burial fluxes, and CaCO 3 and CaSiO 3 weathering feedbacks on the CO 2 and δ 13 C perturbations. The evolution of IRF(CO 2,a ) for experiment SED 500 (no terrestrial weathering feedback but with marine sediments) is comparable to WEA 500 until simulation year 1000 and, it differs thereafter due to varying weathering rates in WEA 500 (Figs. 1a and 2a). There is no addition of alkalinity to the ocean from enhanced terrestrial weathering in SED 500 , and CO 2,a levels out at a higher value than in WEA 500 , with about 8 % of the perturbation remaining constantly airborne in SED 500 . This value is in agreement with previous studies (e.g., Archer et al., 1998). It is somewhat higher than suggested by Eq. (7) as this equation represents a simplification.
CaCO 3 compensation leads, in both WEA 500 and SED 500 , to a further removal of CO 2,a on multimillennial timescales as seen between 1 and 10 kyr in Fig. 2a. After 10 kyr, IRF(CO 2,a ) decreases further as a result of weathering-burial imbalances in the CaCO 3 flux, also referred to as terrestrial neutralization (e.g., Archer et al., 1998;Colbourn et al., 2015). The e-folding timescale for this terrestrial neutralization from SED 500 is about 11 kyr, which is higher than the value provided by Archer et al. (1998) (8.5 kyr) but in the range (8 to 12 kyr) provided by Colbourn et al. (2015).
Comparing δ 13 C atm of SED 500 and WEA 500 , the evolution is initially comparable, but on multimillennial timescales the isotopic perturbation is, somewhat surprisingly and in contrast to the CO 2 perturbation, removed slower in WEA 500 compared to SED 500 (Fig. 2b). In the case of SED 500 the e-folding timescale is on the order of 20 to 30 kyr. This is much shorter than expected from the mean residence time of carbon in the atmosphere-ocean-land-biosphere system of about 90 kyr (carbon inventory: ∼ 40 000 GtC; burial flux: 0.46 GtC yr −1 ). In simulation WEA 500 , with weathering feedbacks enabled, the long-term removal timescale of δ 13 C atm (∼ 70 kyr) is also shorter than the residence time of carbon in the atmosphere-ocean-land-biosphere system. The isotopic perturbation is larger in the surface ocean than in the ocean mean. In turn, the perturbation is removed more efficiently than expected from the residence time of carbon (see Fig. A1). The isotopic signature of the carbon burial flux is mainly determined by the upper-ocean isotopic signature, which is more depleted than the ocean mean for WEA 500 . In the case of SED 500 this gradient is even stronger and lasts longer ( Fig. A1; see Sect. 3.2 for further explanation). These results highlight the importance of spatial patterns within the ocean for the removal of an isotopic perturbation even for these very long timescales. However, the absolute perturbation in δ 13 C a for these small pulse sizes and with δ 13 C = −24 ‰ is rather small after a few thousand years. For example, in the case of SED 500 , δ 13 C a decreases to below 0.05 ‰ after ∼60 kyr. The interannual variability from the LPX-Bern model is of the same magnitude, and model drift also starts to play a role.
In scenario CLO 500 , we exclude the effect of ocean sediment interactions, and we treat the atmosphere, ocean, and land biosphere as a closed system. Considering ocean mixing times of 1-2 kyr, we would expect the CO 2 perturbation for experiment CLO 500 to reach a steady state within 2 kyr. The slow response timescale of peat carbon in the LPX-Bern model, however, leads to a further decline of CO 2,a over multiple millennia (Fig. 2a). About 15 % of the initial perturbation remains constantly airborne in the closed system model setting as expected from Eq. (5). δ 13 C a levels out around 2 kyr at around 3 %, which is as expected from 432 A. Jeltsch-Thömmes and F. Joos: Evolution of atmospheric carbon and δ 13 C perturbations Eq. (6). The further decrease on multimillennial timescales results from model drift in the LPX-Bern model.
Next we compare WEA 500 , where LPX-Bern is used as land model, with 4box 500 , where a four-box land biosphere is used. Differences in IRF(CO 2,a ) between WEA 500 and 4box 500 are due to the different land biosphere models. The four-box land biosphere model takes up less of the CO 2 perturbation than LPX-Bern. The root-mean-square deviation in IRF(CO 2,a ) from the two simulations is less than 1 % of the perturbation. Differences are most pronounced between simulation year 1 and 10 000 (Fig. 2a). Differences in IRF(δ 13 C a ) (Fig. 2b) are largest in the first decades after the pulse, and the root-mean-square deviation in IRF(δ 13 C a ) amounts to less than 1 % between the two simulations.
Finally, we compare the simulations LGM 500 , forced by glacial boundary conditions, and 4box 500 , forced by preindustrial conditions. Differences between the two simulations are solely due to the difference in climate and CO 2 forcing. IRF(CO 2,a ) is lower by a few percent in LGM 500 compared to 4box 500 (Fig. 2a). The difference is largest in the first 10 kyr, and it is up to 12 % of the perturbation. The rootmean-square deviation between the two simulations amounts to ∼ 2 %. The larger oceanic carbon uptake in LGM 500 compared to 4box 500 can be explained by the higher alkalinity simulated for LGM than for preindustrial conditions. A higher alkalinity implies lower values of the Revelle factor ξ and a larger carbon uptake capacity of the ocean. IRF(δ 13 C a ) shows a similar pattern with lower values for the LGM background compared to the preindustrial background (cf. dotted black and dotted dark-violet lines in Fig. 1b). The difference can be understood with the help of Eq. (6); inserting the small glacial value for the atmospheric carbon inventory (381.6 GtC) in the first term yields IRF ∞ (δ 13 C a ) = 0.027 compared to IRF ∞ (δ 13 C a ) = 0.032 when using the preindustrial value (589.4 GtC) in the equation. These estimates of IRF(δ 13 C a ) are in good agreement with the results from LGM 500 and 4box 500 (Fig. 2b). In summary, the influence of glacial compared to preindustrial boundary conditions on IRF(δ 13 C a ) appears modest in our model.

The budgets of the carbon and carbon isotope perturbations
Next, we discuss the budgets of the carbon and 13 C perturbations. Figure 3 shows results for simulation WEA 250 , with a pulse release of 250 GtC, and for simulation WEA −250 , with a pulse removal of 250 GtC. Differences between these two runs are clearly visible in Fig. 3, but the general evolution of the budgets is similar. In the following, we give numerical results for simulation WEA 250 . After 100 years, about 82 GtC is still airborne, and the land and the ocean carbon inventories have increased by about 65 and 105 GtC, respectively (Fig. 3a). As time evolves, the oceanic carbon perturbation grows at the cost of the atmosphere and land. On multimillennial timescales, carbonate compensa- tion and enhanced terrestrial carbonate and silicate weathering increase ocean alkalinity, and thereby the uptake capacity for atmospheric CO 2 and the remaining atmospheric CO 2 fraction decreases. Nevertheless, the carbon perturbation in the atmosphere-ocean-land system exceeds the 250 GtC of the pulse as more carbon is added to the ocean through increased CaCO 3 weathering and redissolution or reduced burial of marine sediments (hatching in Fig. 3a). The maximum perturbation is reached after ∼ 16 kyr with a total of about 381 GtC in the atmosphere-ocean-land system. Thereafter, this inventory decreases due to weathering-burial imbalances. After 100 kyr, about 5 GtC of the perturbation remains airborne, 8 GtC is stored in the land biosphere, and the ocean contains 163 GtC more carbon than before the perturbation. Accordingly, 74 GtC of the initial perturbation is lost from the atmosphere-ocean-land system through increased marine burial compared to terrestrial weathering input (hatching extends below the 250 GtC line in Fig. 3a).
The isotopic perturbation budget is established in units of gigatonnes of carbon per mill (GtC ‰). Any isotopic inven-tory is defined as the product of the corresponding carbon inventory in gigatonnes of carbon (GtC) multiplied by its isotopic signature in ‰. The role of the land biosphere in removing the δ 13 C a perturbation is clearly visible (Fig. 3b). After 100 years, the ocean contains ∼50 % of the isotopic perturbation (in GtC ‰), the land biosphere about 35 %, and the atmosphere less than 15 %. On millennial timescales, the perturbation is slowly removed from the atmosphere-oceanland system. The δ 13 C perturbation of the burial flux (combined POC and CaCO 3 ) follows the negative δ 13 C perturbation in surface DIC. Therefore, the burial flux is what ultimately removes the δ 13 C perturbation from the atmosphereocean-land-biosphere system (Fig. 3b).
We further illustrate the contribution of the burial and weathering fluxes to the carbon and isotope budgets for the simulation with (WEA 500 ) and without weathering feedbacks (SED 500 ) in Fig. 4. The inventory in the relatively fastexchanging pools -atmosphere, ocean, reactive sediments, land biosphere -increases by ∼ 350 GtC in addition to the pulse input of 500 GtC over the course of simulation SED 500 , leading to a total increase of 850 GtC. The increase, in addition to the pulse, is mainly due to a reduction in the burial of CaCO 3 . In contrast, CaCO 3 burial is enhanced in simulation WEA 500 , and ∼ 430 Gt of carbon is removed from the fast-exchanging reservoirs. This removal is only partly compensated by enhanced weathering (∼ 380 GtC), leading to a reduction in the inventory of the fast-exchanging pools from 500 GtC immediately after the pulse to 453 GtC at 100 kyr. Changes in POC burial fluxes are negligible for the carbon budget.
Turning to the isotopic perturbation budget, the initial perturbation due to the pulse input of −12 000 GtC ‰ is mainly removed by anomalies in the POC burial flux both in SED 500 (∼ 10 000 GtC ‰) and WEA 500 (∼ 7000 GtC ‰). Changes in CaCO 3 burial and weathering mitigate about a quarter (∼ 2900 GtC ‰) of the initial isotopic perturbation in the two simulations. This leaves a relatively small isotopic inventory perturbation of about +1100 and −1600 GtC ‰ in the fastexchanging reservoirs at 100 kyr in simulations SED 500 and WEA 500 .
The contributions of the anomalies in the POC and CaCO 3 burial fluxes may be approximately attributed to changes in the signature ( δ) and to changes in flux ( F ) relative to the initial flux, F 0 , and signature, δ 0 , before the perturbation: The isotopic perturbation of the POC burial flux is mainly mediated by a change in the signature of the organic carbon buried ( δ 13 C(POC)). For example, the signature of the total POC burial flux is on average 0.29 ‰ more negative than the corresponding input from weathering in simulation WEA 500 . This change in signature of the total POC burial flux leads to an effective reduction of the isotopic perturbation. The change in POC burial flux ( F (POC)) is, as mentioned above, negligible. Both changes in the signature of the total CaCO 3 burial flux and changes in CaCO 3 burial contribute significantly to the isotopic perturbation. At the end of simulation WEA 500 , the mean δ 13 C signature of CaCO 3 burial is on average 0.14 ‰ more negative than the input from weathering. In turn, the term F 0 · δ for CaCO 3 burial contributes ∼ −3100 GtC ‰ to the mitigation of the initial isotopic perturbation. This is partly counteracted by excess burial of CaCO 3 , with F · δ 0 equal to ∼ 1300 GtC ‰. This yields an overall burial contribution of ∼ −1800 GtC ‰. The δ 13 C perturbation is further removed by increased weathering input of CaCO 3 (see "Model and experimental setup"), contributing about 1100 GtC ‰ over 100 years in WEA 500 . In summary, the isotopic perturbation is mainly mitigated by a change in the signature of the mean POC burial flux. The smaller contributions from the CaCO 3 cycle to the isotopic budget are related to both a change in the signature of the mean CaCO 3 burial flux and to changes in burial and weathering carbon fluxes.

434
A. Jeltsch-Thömmes and F. Joos: Evolution of atmospheric carbon and δ 13 C perturbations

Influence of pulse size
The IRF for atmospheric CO 2 and δ 13 C is sensitive to the magnitude of the pulse size (Fig. 1a, b), and we next address related differences between simulations. For CO 2 , the remaining atmospheric fraction is higher for larger pulse sizes than for smaller pulse sizes. Differences between different small pulses (−250 to 500 GtC) are most significant in the first decades; differences between the large (WEA 5000 ) and any small pulse (Fig. 1a) are most significant in the first millennium. Several processes contribute to these differences.
First, nonlinearities in the carbonate chemistry contribute to differences in IRF(CO 2,a ) in the first centuries to millennia, where the perturbation is largely contained in the atmosphere-ocean-land-biosphere system. We see from Eq. (5) that a controlling factor is the oceanic buffer factor. As pulse size increases, so does the buffer factor, leading to higher values of IRF(CO 2,a ) for WEA 5000 compared to smaller pulse sizes (Fig. 1a).
Second, productivity on land depends nonlinearly on CO 2 (Farquhar et al., 1980), and the loss of carbon from the land biosphere in response to a negative pulse is larger than the uptake in response to a positive pulse of the same size. This contributes to the initially smaller remaining atmospheric fraction for negative pulses (see WEA 250 and WEA −250 in Figs. 1a and 3a). The difference is most strongly pronounced in the first centuries and millennia, but it prevails until the end of the simulation.
Third, the response in marine export production is also nonlinear for pulses of same absolute size but different sign. The initial decrease in export production after a positive pulse is smaller than the increase after a negative pulse of the same size, thereby affecting surface ocean CO 2 concentrations, and thus air-sea gas exchange, and partly leveling out the effect from the land biosphere.
For δ 13 C a , the perturbation is smaller for negative compared to positive pulses of the same magnitude (Fig. 1). The difference can be understood with the help of Eq. (6). Evaluation of Eq. (6, see Appendix) yields an IRF(δ 13 C a ) of ∼ 2 % for the pulse addition of 250 GtC (WEA 250 ) but of only ∼ 1 % for a corresponding removal (WEA −250 ).

Radiocarbon
Finally, we discuss the 14 C perturbation experiments. In experiment 14 C only , a positive pulse of 14 C is added to the atmosphere, but there is no perturbation in CO 2 , analogous to a radiocarbon production pulse by atomic bomb tests or cosmic rays. The initial evolution of the IRF( 14 C a ) is very similar to the one of IRF(δ 13 C a ) (Fig. 1b, c). On longer timescales, radioactive decay additionally removes the perturbation, so that after ∼ 20 kyr less than 0.1 % of the perturbation is left airborne. In simulations 14 C −500 and 14 C dead , a radiocarbon-depleted carbon source of 100 GtC is added to the atmosphere causing a negative perturbation in 14 C and a positive perturbation in CO 2 in the atmosphere. As visible from Fig. 1a for WEA 100 , a small percentage of the CO 2 perturbation remains airborne even after 100 kyr. This explains the persistence of a small atmospheric perturbation in 14 C, well beyond the lifetime of the initially added 14 C. IRF( 14 C a ) is larger for 14 C −500 than for 14 C dead on these very long timescales. This is because the perturbation in CO 2 , and in turn the long-term perturbation in 14 C a , is the same for both experiments. However, the initial perturbation in 14 C -the denominator of the response function (see Eq. 2) -is smaller in 14 C −500 than in 14 C dead .

The δ 13 C perturbation of marine dissolved inorganic carbon and its approximation by PC-EOF
The ocean mean perturbation in δ 13 C of DIC ( δ 13 C DIC ) increases in absolute magnitude during the first millennium, as the atmospheric perturbation is propagated to the ocean (Fig. 5, thick gray line). The most negative mean perturbation (−0.23 ‰ for experiment WEA 500 ) occurs ∼ 1500 years after the pulse release (Fig. 5), reflecting ocean mixing timescales. Thereafter, the δ 13 C DIC perturbation decreases in magnitude as a result of weathering-burial imbalances (see Sect. 3.1, Figs. 3 and 5). At the end of the simulation, less than 15 % of the maximum perturbation (−0.03 ‰) remains in the ocean.
Deep ocean invasion of the δ 13 C DIC perturbation is first accomplished in the North Atlantic Deep Water (NADW) region. After 100 years, the signal has already propagated to a depth of ∼ 3.5 km in the North Atlantic, whereas it is confined to the upper ocean (> 1 km) in the rest of the ocean (Fig. 6a). At year 1000, the δ 13 C DIC perturbation is rather uniformly distributed in the ocean. Only the deep North Pacific (region of oldest waters; ideal age > 1 kyr in Bern3D) and waters at intermediate depths in the Southern Ocean show a smaller perturbation (Fig. 6c). On multimillennial timescales, δ 13 C DIC decreases throughout the ocean. The perturbation, however, is not entirely uniformly distributed as expected from a completed ocean mixing by physical transport. More negative δ 13 C DIC values are present in the upper ocean and in North Atlantic Deep Water compared to the intermediate and deep Southern Ocean, in Antarctic Bottom Water in the Atlantic, and in the deep Pacific (Fig. 6e, g). This may be due to changes in the fractionation during airsea gas exchange at higher temperatures, leading to more negative δ 13 C signatures in the surface ocean. Additionally, gross air-sea fluxes are increased as a result of the higher CO 2 concentrations, decreasing the air-sea disequilibrium. This leads to more negative δ 13 C signatures in mid-and lowlatitude surface waters and less negative signatures in cold high-latitude waters (see Menviel et al., 2015). This surfaceto-deep δ 13 C DIC gradient is dampened in the case of experiment WEA 500 by enhanced weathering of CaCO 3 on land, adding carbon with a positive δ 13 C signature (∼ 2.9 ‰) to the surface ocean. As discussed above, these spatial δ 13 C Figure 5. Modeled global-average perturbation of the isotopic signature of dissolved inorganic carbon in the ocean (δ 13 C DIC ; thick gray line) in response to a carbon pulse to the atmosphere of 500 GtC with a δ 13 C signature of −24 ‰. The thin colored lines show the perturbation as reconstructed with one (blue), two (green), and three (orange) principal components.
gradients are important and co-govern the removal rate of the δ 13 C perturbation. The larger-than-average δ 13 C perturbation in the surface is imparted to phytoplankton and zooplankton, feeding the burial flux of biogenic material. Further, fractionation during phytoplankton growth at elevated CO 2 concentrations is higher, leading to lower δ 13 C values of the POC export flux. This accelerates the removal of the δ 13 C perturbation in comparison to a uniformly mixed ocean.
Next, the temporal evolution of the δ 13 C DIC perturbation and its spatial pattern are described by PC-EOF (Fig. 7). The first PC-EOF pair captures approximately the δ 13 C DIC evolution after 1 kyr, while the second and third PC-EOF pairs capture, together with the first pair, the decadal-tocentennial-scale penetration of the perturbation into the ocean (Fig. 5). The first EOF pattern (Fig. 7a) shows negative values throughout the ocean, with a pattern similar to that modeled for year 50 000 (Fig. 6g). The corresponding PC time series increases during the first decades after the pulse, followed by a slow decline on multimillennial timescales (Fig. 7b). This first PC captures the evolution of ocean mean δ 13 C DIC well after 1 kyr (Fig. 5, blue versus gray lines). In contrast to the first EOF, the second EOF has a clear dipole pattern with negative values in surface and intermediate waters and in the NADW region and positive values below ∼ 1 km depth (Fig. 7c). The corresponding PC time series has high positive values in the first decades of the simulation, declining over the next centuries, becoming negative around year 500, and approaching zero on multimillennial timescales. In combination with the first PC-EOF pair this leads to strongly negative δ 13 C values in surface and intermediate waters in the first centuries of the simulation, whereas on longer timescales the depth gradient in δ 13 C DIC disappears. The third PC-EOF, finally, is of importance during the first century (Fig. 7b), and it captures regional features, e.g., related to North Atlantic Deep Water or intermediate waters (Fig. 7d).
The global-mean δ 13 C DIC is well described by the first two PCs (Fig. 5, green vs. gray line). The third PC-EOF pair mostly adds performance during the first century (cf. green and orange lines in Fig. 5), and it reduces the root-meansquare error (RMSE) between reconstructed and modeled global-mean δ 13 C DIC from 0.003 ‰ with two PC-EOF pairs to 0.001 ‰ with three PC-EOF pairs used in the reconstruction.
For the reconstruction of the spatiotemporal evolution of δ 13 C DIC , we will now use the first three PC-EOF pairs. The RMSE between the modeled and reconstructed 3-D field is highest (0.1 ‰) for the first 3-D output at year 10, after the pulse, and decreases to 0.03 ‰ at year 100. As visible from Fig. 6a and b, largest deviations between modeled and reconstructed δ 13 C DIC are found in the thermocline of the Pacific and in the deep North Atlantic. Further, the modeled perturbation has not yet propagated to the deep Pacific and Atlantic by year 100 (Fig. 6a), whereas the PC-EOF reconstruction displays a δ 13 C DIC signal in these water masses. The performance of the reconstruction increases over the course of the simulation. The RMSE of the δ 13 C DIC fields amounts to 0.01 ‰ at 1 and 10 kyr, and it decreases to 0.007 ‰ at 50 kyr after the pulse. Overall, reconstructing δ 13 C DIC with three PC-EOF pairs shows good results, and deviations from the modeled δ 13 C DIC are of the order of 10 % of the perturbation and smaller on centennialto-millennial timescales (Fig. 6d, f, h). To reconstruct the modeled evolution of δ 13 C DIC precisely during the first few centuries after a perturbation, more than three PC-EOF pairs are necessary. In summary, the spatiotemporal evolution of δ 13 C DIC over ∼ 1000 years to 100 kyr after a pulse input of carbon into the atmosphere is described accurately and conveniently by three PC-EOFs pairs.

Discussion and conclusion
The aim of this study is to investigate the evolution of a δ 13 C perturbation in the Earth system in response to a carbon input to the atmosphere on timescales from centuries to 100 kyr. To this end, the response of the Bern3D-LPX model was probed by pulse-like emissions of carbon to the atmosphere from an external source such as land biosphere carbon release or fossil fuel burning. The δ 13 C response is compared to that of atmospheric CO 2 and radiocarbon.
The impulse response (or Green's) functions (IRF) for atmospheric CO 2 and its δ 13 C signature are fitted by a sum of exponential functions. Additionally, we show that the spatiotemporal evolution of a δ 13 C DIC perturbation is reasonably represented with three PC-EOF pairs. Deviations between the PC-EOF reconstruction and modeled δ 13 C DIC results are largest in the first decades after the pulse, and they are on the order of 10 % of the perturbation on millennial timescales. This allows for the construction of a computationally efficient substitute model, which can be applied to  investigate, for example, reconstructed δ 13 C DIC variations from marine sediments in future studies. The evolution of the atmospheric CO 2 perturbation in the coupled Bern3D-LPX model is comparable to results from earlier studies addressing the pulse-like input of carbon into the atmosphere (e.g., Maier-Reimer and Hasselmann, 1987;Siegenthaler and Joos, 1992;Archer et al., 1998Archer et al., , 2009Ridgwell and Hargreaves, 2007;Colbourn et al., 2015;Lord et al., 2016). The inclusion of peat lands in the LPX-Bern model adds a (small) additional sink on multimillennial timescales, previously not considered in such simulations. For the carbon isotopes, corresponding long-term pulse response simulations that consider ocean-sediment and weathering feedbacks are lacking in the literature. The "short-term" response from years to one or two millennia is consistent with earlier results .
We quantify the processes leading to the different removal timescales of atmospheric CO 2 and δ 13 C atm perturbations with the help of factorial experiments. The role of different processes is further illustrated by analytical descriptions derived from a two-box atmosphere-ocean model and three-box atmosphere-ocean-land model provided in the Appendix. The removal of an atmospheric perturbation in the isotopic ratio 13 C/ 12 C is not buffered by marine carbonate chemistry, unlike CO 2 . The atmospheric δ 13 C perturbation is partitioned among the ocean, land biosphere, and atmosphere roughly in proportion to their carbon inventories over about 2 kyr. This leaves only a very small fraction airborne (1 % to 3 % for pulses of a few hundred gigatonnes of carbon, GtC) after the signal is mixed within these three reservoirs. In contrast, a substantial fraction of carbon emissions and the CO 2 perturbation remain airborne for millennia. On timescales of a few millennia to 20 kyr, the initial atmospheric CO 2 perturbation is lowered by carbonate compensation to a fraction of about 8 %. In contrast, this chemical buffering hardly affects the atmospheric δ 13 C perturbation. On even longer timescales, the CO 2 and the isotopic perturbation are removed by different, though related, processes. For δ 13 C, the continuous flux of biogenic calcium carbonate and organic carbon particles carries the isotopic perturbation from the surface ocean to the lithosphere. For CO 2 , the excess weathering of silicate rocks in response to the CO 2 and climate perturbation adds alkalinity to the ocean, which leads to a complete removal of the atmospheric CO 2 perturbation.
Gradients in the δ 13 C DIC perturbation strongly influence the long-term removal timescale of the isotopic perturbation. In the pulse experiments, δ 13 C is enlarged in the surface relative to the deep ocean by temperature-dependent fractionation during air-sea gas exchange (Mook, 1986) and altered air-sea disequilibrium resulting from changes in gross exchange fluxes due to altered CO 2 . Changes in CaCO 3 weathering, on the other hand, diminish the surface perturbation (e.g., experiment WEA 500 vs. SED 500 in Fig. 1a). All these processes lead to spatial gradients in the δ 13 C perturbation in 438 A. Jeltsch-Thömmes and F. Joos: Evolution of atmospheric carbon and δ 13 C perturbations the ocean. As a result of these gradients in δ 13 C, the isotopic perturbation is removed faster by particle burial than expected from a uniformly mixed ocean. This highlights the importance of resolving spatial structures in the perturbation to represent its removal timescales. More generally and beyond pulse-release experiments, any process that changes the δ 13 C signature of the relentless flux of organic and calcium carbonate particles will influence the δ 13 C loss flux from the ocean-atmosphere-land-biosphere system to the lithosphere (for examples, see Jeltsch-Thömmes et al., 2019;Roth et al., 2014;Tschumi et al., 2011).
Our results have consequences for the interpretation of the difference in δ 13 C between similar climate states such as the Penultimate Glacial Maximum (PGM) and the Last Glacial Maximum (LGM). Substantial temporal δ 13 C differences between these periods are recorded in ice cores (Schneider et al., 2013;Eggleston et al., 2016) and in marine sediments (Hoogakker et al., 2006;Oliver et al., 2010). Different mechanisms -such as changes in the δ 13 C signature of weathering and burial fluxes, varying contribution of volcanic outgassing of CO 2 , and changes in the amount of carbon stored in the land biosphere, especially in yedoma and permafrost soilshave been discussed for the PGM-LGM δ 13 C offset (e.g., Lourantou et al., 2010;Schneider et al., 2013). An internal reorganization of the marine carbon cycle without considering changes in burial may not explain the offset. With no changes in the weathering-burial balance, the mass in the atmosphere-ocean-land system remains constant. Then, a change in atmospheric δ 13 C would require an opposing δ 13 C change in the ocean. This appears to be in conflict with marine and ice core records, which suggest that δ 13 C increased both in the ocean and in atmosphere by about 0.3 ‰ between the PGM and LGM (Hoogakker et al., 2006;Oliver et al., 2010;Schneider et al., 2013;Eggleston et al., 2016). However, internal marine carbon-cycle reorganizations, e.g., due to changes in circulation, may have altered δ 13 C in the surface ocean, and in turn δ 13 C of the burial flux, and thereby the balance between weathering input and burial.
The Bern3D results suggest that changes in organic carbon storage on land (or in the ocean) are not a likely explanation of the isotopic offset between the PGM and LGM. The δ 13 C signal of a transfer of organic carbon of plausible magnitude to the atmosphere and ocean would have been attenuated too much over time to cause the observed PGM-LGM δ 13 C difference. Schneider et al. (2013) estimated required changes in land carbon storage to match the ice core and marine temporal offsets in δ 13 C using preliminary Bern3D results for pulse-release simulations. Our results show, in agreement with these earlier results, that the isotopic perturbation associated with a terrestrial carbon release is attenuated to less than 15 % within two millennia for pulse sizes of up to 5000 GtC and declines thereafter (see black line in Fig. 1b). This would require several thousand GtC to have been stored additionally in the land biosphere during the LGM compared to the PGM in order to explain the δ 13 C off-set (see also Schneider et al., 2013). This amount seems large in light of estimated total carbon stocks of ∼ 1500 GtC in perennially frozen soils in Northern Hemisphere permafrost regions today (Tarnocai et al., 2009) and an even smaller inventory at the LGM (Lindgren et al., 2018).
Taken together, these ad hoc considerations suggest that long-term imbalances in the weathering (including volcanic emissions) and burial fluxes appear to be a plausible cause for the δ 13 C differences between PGM and LGM. However, further work, which considers the transient evolution of CO 2 , δ 13 C, and other proxies, is required to gain further insight into the contributions of individual mechanisms to long-term δ 13 C changes recorded in ice and marine cores.
Further, the different behavior of CO 2 and δ 13 C atm perturbations also has consequences for centennial-scale CO 2 and δ 13 C atm variations such as during Heinrich stadials 4 and 1 (HS4 and HS1). Variations in CO 2 and δ 13 C atm during HS4 and HS1 have been measured on Antarctic ice cores. The data show an increase in CO 2 on the order of ∼ 10 ppm over 200-300 years, accompanied by a decrease in δ 13 C atm of ∼ 0.2 ‰. This results in a value of roughly −0.02 ‰ ppm −1 for the ratio of the change in δ 13 C to the change in CO 2 (r = δ 13 C/ CO 2 ). A release of terrestrial carbon has been discussed as a possible cause of these events (Bauska et al., 2016(Bauska et al., , 2018.
We utilize our pulse response simulations under PI and LGM boundary conditions to estimate the changes in r in response to a transient terrestrial carbon input with an isotopic signature of −24 ‰. The results for r will depend on the evolution of the emissions into the atmosphere. For a pulselike input at time zero, the evolution of r after the pulse is shown in Fig. 8a for the first 2000 years. r is about −0.08 (−0.06) ‰ ppm −1 immediately after the release under LGM (PI) boundary conditions. The difference between LGM and PI values of r comes from the difference in the atmospheric C and 13 C inventory to which the pulse is added. The magnitude of r is reduced rapidly within the first 25 years and remains approximately constant at −0.015 (−0.01) ‰ ppm −1 thereafter. Thus, r recorded in ice after a hypothetical pulse release of carbon would be above −0.02 ‰ ppm −1 if the age distribution of the ice core samples is wider than a few decades. However, emissions may occur more gradually. For simplicity, we assume an emission pulse sustained over period T according to a two-sided Heaviside function; terrestrial emissions, e(t), suddenly increase to unity, remain elevated for a period T , and are then immediately reduced again to zero. r is then readily calculated from the IRFs (r = T 0 dt · e(t)· IRF(δ 13 C)(T − t)/ T 0 dt · e(t)· IRF(CO 2 )(T − t)). The results for r are again plotted in Fig. 8a. r is decreasing with the duration T of sustained emissions and falls below −0.02 ‰ ppm −1 for T larger than about 90 (30) years for LGM (PI) boundary conditions. Finally, we also performed simulations with the Bern3D-LPX under PI boundary conditions to account for nonlinearity that may not be captured in the IRF representation. The releases of varying amounts of Carbon is released as a single pulse at year 0 (green lines) or assumed to be released at a uniform rate from year 0 to the year plotted (purple lines; see main text). The dashed horizontal line indicates a δ 13 C-to-CO 2 ratio of −0.2 ‰ to 10 ppm as seen in ice core records for HS1 and HS4. (b) Response in CO 2 and δ 13 C atm to a uniform release of 20, 40, 80, 120, and 160 GtC (colors) of isotopically light carbon (δ 13 C = −24 ‰) to the atmosphere over 100, 200, 300, and 400 years (numbers). Shown is the mean over 30 years starting in the last year of the release. The red cross indicates the approx. δ 13 C-to-CO 2 change as seen in ice core records for HS1 and HS4 (Bauska et al., 2016(Bauska et al., , 2018. terrestrial carbon with δ 13 C = −24 ‰ over 100 to 400 years yield a linear response in CO 2 and δ 13 C atm (Fig. 8b) with a slope r of about −0.01 ‰ ppm −1 , in agreement with the IRF results.
There are uncertainties in these estimates of r. We assumed that the emitted carbon was about −18 ‰ depleted with respect to the atmosphere. The relative uncertainty in this value translates directly to the relative uncertainty in r. The release of material that is isotopically even more depleted would bring our model results into better agreement with the ice core estimate of r. Observations and land model simulations suggest that the δ 13 C depletion in vegetation and soils may vary between −24 ‰ for heavily discriminated C 3 plant material and to −5 ‰, which is indicative of C 4 plant material (Keller et al., 2017). C 4 plants were likely more abundant during glacials than interglacials (Collatz et al., 1998). Also, the values shown in Fig. 8b for century-scale emissions are 30-year means. A wider gas age distribution would reduce the perturbation in δ 13 C atm more than in CO 2 (see Köhler et al., 2011). On the other hand, the CO 2 increase may have occurred in less than 100 years (Bauska et al., 2018); correspondingly, a more negative value of r is expected for such a fast (< 100 years) increase. Given these uncertainties in our model estimates, as well as in the determination of r from the ice core measurements, we are not in a position to firmly conclude on whether a terrestrial carbon release was the sole mechanism for the HS1 and HS4 CO 2 events or not. Other mechanisms in addition to a terrestrial release may have contributed to the CO 2 and δ 13 C atm excursions during HS4 and HS1.
We acknowledge that most of the idealized simulations presented in this study were run under preindustrial boundary conditions. A limited set of sensitivity simulations suggests that the IRFs for CO 2 and δ 13 C are similar under preindustrial and glacial CO 2 in our model. The direct influence of lower atmospheric CO 2 on IRF(δ 13 C) appears limited. Further work is needed to better quantify the dependency of IRFs on the climate state.
In conclusion, the responses of atmospheric CO 2 and of δ 13 C atm to a carbon input (or removal) are different, with implications for the interpretation of paleoproxy data. The difference in response is related to the aquatic carbonate chemistry and, on millennial-to-multimillennial timescales, to the burial weathering cycle, which affect CO 2,a and its isotopic ratio differently. Our study highlights the importance of carbon isotope fluxes to and from the lithosphere on timescales when the ocean-atmosphere-land-biosphere system has so far often been assumed to represent a closed system. Perturbations in δ 13 C are removed from or modified within the atmosphere-ocean-land-biosphere system by the continuous flux of biogenic organic and calcium carbonate particles buried in marine sediments. This isotopic burial flux is substantially influenced by spatial gradients in δ 13 C within the ocean, calling for a spatially resolved representation of the marine carbon cycle to address isotopic change. Carbon isotope fluxes to and from the lithosphere need to be considered for the interpretation of changes in atmospheric and marine δ 13 C during glacial-interglacial, and longer, periods. The difference in the IRF for carbon and the isotopic ratios can be understood in terms of the aquatic carbonate chemistry. We develop approximate expressions for the IRF of atmospheric CO 2 (Eq. A7) and δ 13 C (Eq. A29), starting in this section with IRF(CO 2 ).
For small pulses, changes in sea surface temperature, in CO 2 solubility, in ocean circulation, and in the marine biological cycle are small and neglected in the following. Changes in land biosphere carbon storage can be considerable, in particular during the first century. However, the role of changes in land carbon stock becomes smaller when more carbon is taken up by the ocean, and thus it is neglected for simplicity.
We consider a two-box model of the atmosphere (index: a) and ocean (index: o) with area A = A a = A o and heights h a and h o . All boxes are assumed to be well mixed. Before the pulse addition, the model is in equilibrium and the initial CO 2 partial pressure is pCO 2,0 =pCO 2,a,0 =pCO 2,o,0 . The initial concentration of dissolved inorganic carbon (DIC) in the ocean is DIC 0 , and the carbon inventory in the ocean is given by N o = DIC · h o · A. We relate the atmospheric car-bon inventory, N a , to a time-varying atmospheric concentration C a and to the fixed volume of the model atmosphere: N a = C a · h a · A. We are free to select the units of the modelspecific concentration C a . Here, we express C a in units of DIC 0 and set the initial atmospheric concentration C a,0 equal to DIC 0 . The atmospheric inventory is proportional to the atmospheric partial pressure, pCO 2,a , and, correspondingly, C a has to be proportional to pCO 2,a . This leads to the following definition of C a : where C a is specifically defined for our box model and should not be confused with the concentration or mixing ratio of CO 2 in the real-world atmosphere. The atmospheric scale-height h a is given by the requirement that volume times concentration equals inventory: The perturbation ( ) from the initial equilibrium in pCO 2 , and similarly for other quantities, is defined by pCO 2 (t) = pCO 2 (t) − pCO 2 (t 0 ). The perturbation in pCO 2,o is related to that in DIC by the Revelle factor ξ (Revelle and Suess, 1957): We note that the carbonate chemistry is nonlinear, and ξ varies with environmental conditions (DIC, alkalinity, temperature, salinity). Assuming a constant Revelle factor, ξ is therefore an approximation. In particular, ξ increases with increasing DIC. As a consequence, the "effective" value of ξ increases with increasing carbon emission, and thus it is larger when the magnitude of the carbon pulses is larger. This yields, for the perturbation in the net air-to-sea carbon flux per unit area, f a→o,net : where pCO 2,a and pCO 2,o were replaced with the help of Eqs. (A1) and (A3) to obtain the right-hand side of Eq. (A4). k g is the gas transfer rate relative to the partial pressure (Wanninkhof, 1992), and g is the rate relative to DIC: The IRF, or equivalent airborne fraction, is defined by for t > t 0 .
P is the magnitude of the initial carbon pulse at time t 0 , and this released carbon is conserved within the model atmosphere and ocean. The net air-to-sea flux becomes zero when a new equilibrium at time t ∞ is reached; it follows from Eq. (A1) that C a,0 = DIC 0 and from Eq. (A4) that DIC ∞ = C a,∞ /ξ . Finally, using this, as well as N a = C a · h a · A and N o = DIC·h o ·A = C a /ξ ·h o ·A, yields the following for the IRF in an ocean-atmosphere system at equilibrium: Distribution of a δ 13 C perturbation within the atmosphere-ocean-land-biosphere system We now turn to the isotope 13 C. The isotopic ratio is given by the following (see Mook, 1986): 13 R(C)= 13 C/ 12 C.
And the isotopic signature in per mill is given by 13 R SD is the standard ratio. Isotopic discrimination is described by factors using the symbol α or by additive terms in δ notation using symbol ε in per mill. It holds true that In the following, we neglect for simplicity the approximately 1 % difference between the concentration of 12 C and total C. (This difference may be accounted for in the definition and values of the isotopic discrimination factors.) We address first the response during the first two millennia after the pulse. In the following paragraphs, we will develop a simple analytical expression for the IRF (δ 13 C a ) after a new equilibrium for the atmosphere-ocean-land-biosphere system is reached. Thus, we neglect carbonate compensation and weathering and burial fluxes. We consider exchange between the atmosphere and ocean and between the atmosphere and land biosphere reservoirs, and we assume that these three reservoirs are well mixed.

The box model
The carbon and 13 C budgets for the three-box model are described by where N is the reservoir size (e.g., in GtC) and indices a, o, and b denote the atmosphere, ocean, and land biosphere boxes, respectively. Index i refers either to carbon or 13 C. F denotes fluxes between reservoirs. δ is the Kronecker symbol, and P i is the pulse released to the atmosphere at time t = t 0 = 0.
The initial atmospheric δ 13 C perturbation The IRF(t) for δ 13 C a is given by the perturbation in the isotopic signature, δ 13 C a (t), divided by the initial (ini) perturbation, δ 13 C ini a ; δ 13 C ini a is the perturbation immediately after a carbon input of amount P and with signature 13 R P (δ 13 C P ). Mass balance implies that the amount of 13 C before and after the pulse release is equal: We convert the ratios to δ units and subtract on both sides (N a,0 + P ) · δ 13 C a,0 to get the initial perturbation: δ 13 C ini a = P N a,0 + P · (δ 13 C P − δ 13 C a,0 ).
δ 13 C ini a is proportional to the difference between the signature of the pulse and the initial atmospheric signature. Thus, the initial perturbation is relative to the atmospheric δ 13 C signature.
The relationship between δ 13 C a , δ 13 C o , and δ 13 C b at equilibrium In this subsection, relationships between the isotopic signatures are developed by considering equilibrium between the three boxes. This will allow us to simplify the isotopic mass balance equations further below.

(A16)
This equation gives the equilibrium relationship between the atmospheric and oceanic signature. α a,o denotes the equilibrium discrimination factor between pCO 2,a and DIC. α a,o depends on temperature and somewhat on carbonate chemistry (Mook, 1986), and it is about 1.008 and thus close to 1. A similar relationship is readily developed for the isotopic ratio of the total carbon in the land biosphere. The gross fluxes between the land and atmosphere are 13 F a→b = F a→b · 13 α a→b · 13 R pCO 2,a 13 F b→a = F b→a · 13 α b→a · 13 R(N b ). (A17) This yields at equilibrium 13 R eq (pCO 2,a ) = 1 α a,b · 13 R eq (N b ).
δ 13 C b,eq = 1 + ε a,b 1000 δ 13 C a,eq + ε a,b ∼ = δ 13 C a,eq + ε a,b . (A19) Equations (A16), (A18), and (A19) allow us to express the isotopic signature of the land biosphere and the ocean by the isotopic signature of the atmosphere at equilibrium. These equations hold true both for the initial equilibrium (t 0 ) before the pulse release and when a new equilibrium is reached at t ∞ . In practice t ∞ is about 2 kyr. It follows immediately from Eq. (A19) that the isotopic signatures at the new equilibrium are about equal for the three reservoirs: In other words, the isotopic signal is mixed uniformly through the atmosphere, ocean, and land biosphere when assuming that these three reservoirs are well mixed as in this appendix section.
Solving the mass balance equations of the isotopic perturbation for ∆δ 13 C The mass balance equations for the ocean-land-biosphereatmosphere system are given for carbon by The second equation describes the perturbation only. We use Eqs. (A16) and (A18) to replace N o · 13 R o with α a,o ·N o · 13 R a , similarly for the land reservoir, and write the mass balance for 13 C.
Now we convert to δ units: We subtract on both side the mass balance equation for carbon (Eq. A21) multiplied by a factor of 1000. Note that α a,o and α a,b are present in the carbon terms in Eq. (A23) but not in the mass balance for carbon (Eq. A21). This difference results in the first two terms in Eq. (A24) below.
Next, we subtract on both sides the third left-hand side term, (N a,0 + α a,o · N o,0 + α a,b · N b,0 ) · δ 13 C a,0 , as well as the term P · δ 13 C a,0 . We rearrange and use again the mass balance for carbon to get − N o,∞ · (δ 13 C DIC,0 − δ 13 C a,0 · (1 + ε a,o /1000)) − N b,∞ · (δ 13 C b,0 − δ 13 C a,0 · (1 + ε a,b /1000)) − P · (δ 13 C a − δ 13 C P ) The terms with ε/1000 represent a small correction (< 0.02) and can be safely neglected. We replace P with N a,∞ + N o,∞ + N b,∞ in the above equation and exchange the left and right side to get the following mass balance for the δ 13 C perturbation: (N a,0 + N o,0 + N b,0 + P ) · δ 13 C a,∞ = N a,∞ · δ 13 C P − δ 13 C a,0 + N o,∞ · δ 13 C P − δ 13 C DIC,0 with Equation (A26) represents the mass balance for the isotopic perturbation and its interpretation is as follows. The term on the left-hand side equals the carbon mass in the system times the perturbation in the signature, and it therefore equals the total isotopic perturbation (GtC ‰). In equilibrium, the δ 13 C perturbation, δ 13 C ∞ , is well mixed (Eq., A20) within the total ocean-atmosphere carbon reservoir (N a,0 + N o,0 + N b,0 + P ). The terms after the first equal sign provide a different view of the total isotopic perturbation. In this picture, a fraction N l,∞ of the pulse is added to each reservoir (index l) with the perturbation in signature given by (δ 13 C P − δ 13 C l ). This is then summarized after the second equal sign by the product of the carbon input P times the isotopic perturbation of the pulse relative to the mean signature of the reservoirs absorbing this additional carbon (δ 13 C P − δ 13 C mean ). By rearranging Eq. (A26) we get the perturbation in the isotopic signature of the atmosphere-ocean-land system: In other words, the perturbation in the signature scales with the carbon mass ratio of the pulse input to the total mass in the system. This perturbation is also proportional to the difference between the signature of the pulse and the "mean" signature of the system as defined by Eq. (A27).
The fraction of the initial perturbation in the atmosphere that remains airborne corresponds roughly to the ratio of the carbon in the atmosphere immediately after the pulse release to the total carbon in the system. The second fraction on the right-hand side modifies this ratio by about +20 % for a typical negative isotopic pulse perturbation applied in this study, because the mean signature is slightly heavier than the atmospheric signature. The numerator of this modifier reflects that the pulse perturbation is distributed within the system at the new equilibrium, while the denominator reflects that the pulse is initially added to the atmosphere.

Ocean invasion: numerical examples
Equations (A7) and (A29) provide simple expressions to approximate the IRF for the perturbation in atmospheric CO 2 and its isotopic signature after a pulse-like carbon input. These equations illustrate the fundamental difference in the IRF for CO 2,a and for δ 13 C a . To ease interpretation, we neglect the land biosphere for the moment. Then, in the limit of P → 0, Eqs. (A7) and Eq. (A29) are formally identical for carbon and 13 C (index: i): However, for δ 13 C a , a Revelle factor of 1 instead of about 12 applies, while the modifier M is exactly 1 for CO 2,a and around 1 for δ 13 C a . The perturbation in the isotopic signature is diluted by the entire carbon inventory in the system, whereas the "available" ocean inventory is reduced by the Revelle factor for the dilution of an atmospheric CO 2 perturbation.
Numerically, the preindustrial carbon inventory is 600 GtC in the atmosphere and 37 400 GtC in the ocean. These values are equivalent to a scale height of 69 m for the atmosphere and of 4300 m for the ocean, when assuming a preindustrial (surface) concentration of DIC of 24 g m −3 and an ocean area of 3.6 × 10 14 m 2 . The buffer factor is about 12 for small perturbations. This yields for IRF ∞ (CO 2,a ) a value of about 0.16. In other words, a fraction of about 16 % of a (small) carbon input into the atmosphere is still airborne after about 1 to 2 kyr, when the ocean and atmosphere have approached a new equilibrium. The observed penetration of chlorofluorocarbons and bomb-produced radiocarbon suggests that an atmospheric perturbation penetrates about the top 300 m of the ocean within a decade. Thus, we expect that about 70 % 444 A. Jeltsch-Thömmes and F. Joos: Evolution of atmospheric carbon and δ 13 C perturbations of the initial perturbation is still found in the atmosphere after a decade. These numbers are comparable in magnitude to the IRF of CO 2 shown in Fig. 1a for 10 years and 2 kyr after the pulse input. The buffer factor increases with the magnitude of the perturbation in DIC and pCO 2 and thus with pulse size P . Hence, IRF increases with increasing pulse size, as again shown in Fig. 1a. We note that any carbon uptake or release by the land biosphere is not taken into account in Eq. (A7) or in this discussion, in contrast to the results discussed in the main text.
To determine the IRF for δ 13 C a , we take into account that the isotopic signal is also entering the land biosphere, and we assume a total inventory in the atmosphere-ocean-land system of 40 000 GtC. This yields an IRF ∞ of 1 % to 3.2 % for P varying between −250 and +500 GtC and of about 18 % for P equal to 5000 GtC and using δ 13 C mean = 4 ‰ (calculated with Eq. A27). These values are in agreement with the estimates shown in Fig. 1b for year 2000. Assuming that the isotopic perturbation has mixed within a layer of 300 m in the ocean (∼ 2600 GtC) and with the living vegetation on land (500 GtC) yields a "decadal" IRF of about 12 %-31 % for small pulse sizes (−250 to 500 GtC), again in agreement with the Bern3D-LPX results.
Carbonate compensation and terrestrial neutralization of the CO 2 perturbation Next, we address CaCO 3 compensation of the carbon added by the pulse. We derive an expression for the IRF(CO 2,a ) at the time when CaCO 3 compensation of the pulse is completed -about 10-20 kyr after the pulse input. The following calculations are based on Archer et al. (1998). As above, we assume the ocean to consist of a single well-mixed box which is, after CaCO 3 compensation, again in equilibrium with the atmosphere. Concentrations of CO 2 , HCO − 3 , and CO 2− 3 in seawater are related by K 1 and K 2 are the apparent dissociation constants for carbonic acid and K H is the Henry's law solubility product for CO 2 . For simplicity, we assume the ratio of K 1 , K 2 , and K H , which are temperature and salinity dependent, to be constant over the course of the experiments. With the above assumption and by taking the ratio of initial and final states, we rearrange to get refers here to the difference between the final and initial state. With pulse sizes substantially smaller than the oceanic DIC inventory, N o N o,0 , we approximate From Fig. 3a we see that on multimillennial timescales the majority of the pulse is taken up by the ocean. Further, CaCO 3 compensation and weathering-burial imbalances in the CaCO 3 cycle (terrestrial neutralization) add (or remove) additional carbon to (from) the ocean. We thus assume that N o ≈ 1.5 · P ; P is the pulse size in GtC. The factor 1.5 is model dependent and slightly different here from the value provided by Archer et al. (1998) or in previous versions of the Bern3D model. This yields for the relative increase in pCO 2 pCO 2 pCO 2,0 ≈ 3 · P N o,0 .
The remaining atmospheric fraction is thus independent of the pulse size for P N o,0 . Its magnitude results directly from the carbonate chemistry (Archer et al., 1998). With a preindustrial ocean DIC inventory of 37 400 GtC and surface ocean pCO 2 of 284 µatm, we get a remaining airborne perturbation of ∼5 %. This is in agreement with results from the Bern3D-LPX model after about 10-20 kyr (Fig. 1a).
Carbonate compensation of the δ 13 C perturbation Carbonate compensation has almost a negligible influence on the δ 13 C perturbation. The carbonate compensation process acts, aside from small isotopic discrimination, equally on the different isotopes. Therefore, 12 CO 2,a and 13 CO 2,a both decrease according to Eq. (A38), leaving the isotopic ratio, to first-order approximation, unchanged. Additional effects are, here formulated for a positive pulse, as follows. The dissolution of CaCO 3 adds an isotopic perturbation to the oceanatmosphere-land system. The amount of carbon added corresponds to about half of the pulse size (P ). CaCO 3 has a δ 13 C signature of about +2.9 ‰ compared to −24 ‰ of the pulse. We reference the CaCO 3 isotopic perturbation relative to the mean signature of the ocean in the model (∼ 0.8 ‰) as this additional carbon remains mainly in the ocean, while the mean signature that applies for the pulse is, according to Eq. (A27), about −3 ‰ for pulse sizes between −250 GtC and 500 GtC. Thus the isotopic perturbation by CaCO 3 dissolution is about 20 times smaller than that of the pulse (0.5×P ×(2.9 ‰−0.8 ‰) : P ×(−24 ‰+3 ‰) = −1 : 20). The dissolution of CaCO 3 enlarges the carbon reservoir by which the initial atmospheric δ 13 C is diluted. However, the amount of added carbon is, except for the 5000 GtC pulse, small compared the total inventory of the system of about 40 000 GtC. In summary, the influence of CaCO 3 compensation on δ 13 C a is small, and Eq. (A29) is still approximately valid for the timescale of CaCO 3 compensation.
Removal of the remaining CO 2 perturbation by silicate weathering Finally, on timescales of hundreds of thousands of years, imbalances in the weathering and burial fluxes remove any remaining perturbation in carbon from the atmosphere-oceanland system. Silicate rock weathering leads to a net removal of carbon from the atmosphere by the following simplified reaction (see Colbourn et al., 2013, for more details): 2CO 2(aq) + H 2 O (l) + CaSiO 3,(s) → Ca 2+ (aq) + 2HCO − 3(aq) + SiO 2(aq) .
The timescale for this silicate weathering feedback has been determined to be on the order of 240 kyr (Colbourn et al., 2015) or 270 kyr (Lord et al., 2016) but with a large spread. Silicate rock weathering removes the remaining atmospheric perturbation (∼ 8 %). Further, the additional alkalinity added to the ocean deepens the saturation horizon, resulting in higher carbon burial in marine sediments, and thus it also removes the DIC perturbation (see Fig. 3a).
Removal of the remaining δ 13 C perturbation by burial of organic matter and CaCO 3 In the case of δ 13 C, silicate rock weathering has no direct effect on the δ 13 C signature. The perturbation is removed on these long timescales by the replacement of carbon through burial of organic matter and CaCO 3 and weathering inputs: d dt (N S · δ 13 C S ) = −F burial · δ 13 C burial , with N S and δ 13 C S equal to the total mass and mean isotopic perturbation in the atmosphere-ocean-land system; F burial being the burial flux of carbon leaving the system; and δ 13 C burial being its isotopic signature. In simulations with constant weathering rates, the residence time of carbon in the atmosphere-ocean-land system is on the order of τ = N/F burial (∼ 90 kyr) (burial or weathering input is 0.46 GtC yr −1 ; total amount: ∼ 40 000 GtC). Surprisingly, in the factorial simulation with sediments and constant weathering fluxes (SED 500 ), the perturbation is removed much faster than expected from the mean residence time (Fig. 1b, dotted  line), because the isotopic perturbation of the burial flux is larger than that of the ocean-atmosphere-land system. The long-term removal rate of the δ 13 C perturbation is sensitive to δ 13 C gradients in the ocean. The isotopic signature of the burial flux is mainly determined by the upperocean signature. Changes in the gross exchange fluxes of carbon between the atmosphere and the ocean in response to perturbed CO 2 and changes in the isotopic fractionation of air-sea fluxes in response to perturbed temperatures increase the δ 13 C perturbation in the upper ocean relative to the deep (Fig. A1). In the case of enabled weathering feedbacks, the additional flux of isotopically heavy carbon from excess weathering partly mitigates surface-to-deep δ 13 C gradients. In turn, the perturbation is removed faster in simulations with constant weathering compared to simulations with enabled weathering feedbacks. Taken together, this highlights the role of spatial patterns of the δ 13 C perturbation even on these very long timescales.

Parameterization of weathering feedbacks on land
We use the 0-D version of the Rock Geochemical Model v0.9 (Colbourn et al., 2013). Weathering of CaCO 3 and CaSiO 3 is parameterized as functions of temperature (T ), runoff (R), and productivity (P ) in the following form for CaCO 3 : and for CaSiO 3 Subscript "0" indicates the initial value diagnosed prior to the pulse; k Ca stems from a correlation between temperature and HCO − 3 ion concentration of groundwater with a value of 0.049; E a is the activation energy for silicate weathering (63 kJ mol −1 ); β is taken as 0.65; and R indicates runoff and is parameterized as a function of temperature: with k run = 0.025. In the case of experiments 4box 500 and LGM 500 , LPX-Bern was not coupled and instead productivity was parameterized, again following equations given in Colbourn et al. (2013): with C being atmospheric CO 2 in parts per million (ppm). Reference values (0) for temperature, runoff, and productivity are diagnosed in the 100 years prior to the pulse. For further details and discussion of parameter value choices, the reader is referred to Colbourn et al. (2013).