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₂ and in its isotopic ratios (δ¹³C, Δ¹⁴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₂ and δ¹³CO₂ 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₂ as explained by aquatic carbonate chemistry. On multimillennial timescales, the CO₂ perturbation is removed by carbonate compensation and silicate rock weathering. In contrast, the δ¹³C perturbation is removed by the relentless flux of organic and calcium carbonate particles buried in sediments. The associated removal rate is significantly modified by spatial δ¹³C gradients within the ocean, influencing the isotopic perturbation of the burial flux. Space-time variations in ocean δ¹³C perturbations are captured by principal components and empirical orthogonal functions. Analytical impulse response functions for atmospheric CO₂ and δ¹³CO₂ are provided.

Results suggest that changes in terrestrial carbon storage were not the sole cause for the abrupt, centennial-scale CO₂ and δ¹³CO₂ variations recorded in ice during Heinrich stadials HS1 and HS4, though model and data uncertainties prevent a firm conclusion. The δ¹³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 fluxes. Our study highlights the importance of isotopic fluxes connected to weathering–sedimentation imbalances, which so far have been often neglected on glacial–interglacial timescales.

1 Introduction

Research efforts to understand the past and future evolution of atmospheric CO₂ 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₂ (δ¹³C_atm) and of dissolved inorganic carbon (δ¹³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 δ¹³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₂ and δ¹³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., 2010, 2017; Schmitt et al., 2012; Bauska et al., 2018). Reconstructed variations in CO₂ and δ¹³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 (e.g., Bauska et al., 2016, 2018); however, other scenarios might be possible.

The perturbations in CO₂ and δ¹³CO₂ 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₂, over the first decades and centuries, this removal is because of ocean and land uptake; on millennial timescales, it is CaCO₃ compensation; finally, on timescales of hundreds of thousands of years, the remaining atmospheric CO₂ perturbation is removed by enhanced silicate rock weathering (e.g., Archer et al., 1998; Joos et al., 2001, 2013; Colbourn et al., 2015; Lord et al., 2016). In the case of δ¹³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₂ 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 reduced-form 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₂ 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₂ concentrations (Hooss et al., 2001). The perturbation in the atmosphere or ocean is then given by the convolution integral of the emission history and the corresponding IRF (e.g., Joos and Bruno, 1996; Joos et al., 1996). 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₂ perturbation. IRFs have been widely used to calculate remaining atmospheric fractions of anthropogenic CO₂ 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; Joos and Bruno, 1996; Archer et al., 1997, 1998, 2009; Ridgwell 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 isotopes of CO₂ are scarce (Schimel et al., 1996; Joos and Bruno, 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₂ 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 (Hooss et al., 2001; 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₂ and its δ¹³C signature in response to an external carbon input. To this end, we performed pulse-release (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 (Δ¹⁴C) are briefly addressed. Implications of our results for the PGM–LGM δ¹³C difference and for past centennial CO₂ variations are discussed. Simple expressions that illuminate the fundamental differences in the removal timescales for an atmospheric perturbation in CO₂ versus the associated perturbation in δ¹³CO₂ 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 ¹⁴C production by cosmic radiation or atomic bomb tests.

2 Model and experimental setup

2.1 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., 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, 1999; Orr et al., 1999, 2017; Wanninkhof, 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).

Marine productivity is calculated as a function of light availability, temperature, and nutrient concentrations (P, Fe, Si), and it is restricted to the euphotic zone in the uppermost 75 m. Tracers, including dissolved inorganic carbon (DIC) and semilabile dissolved organic carbon (DOC); the corresponding isotopic forms; and alkalinity (Alk), phosphate (PO₄), oxygen (O₂), iron (Fe), silicon (Si), and an ideal age tracer are transported by advection, diffusion, and convection. Biogeochemical cycling is described in detail in Tschumi et al. (2011), Parekh et al. (2008), and subsequent studies (e.g., Menviel and Joos, 2012; Menviel et al., 2012; Roth and Joos, 2012; Roth et al., 2014; Menviel et al., 2015; Battaglia et al., 2016; Battaglia and Joos, 2018a, b; Jeltsch-Thömmes et al., 2019). Isotopic fractionation between model components is documented in Jeltsch-Thömmes et al. (2019).

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₃, opal, particulate organic carbon (POC), clay) and seven tracers in the porewater (DIC, DI¹³C, DI¹⁴C, alkalinity, phosphate, oxygen, and silicic acid) are modeled. The dissolution of CaCO₃ and POC depends on the respective weight fraction of CaCO₃ and POC in the solid phase of the sediment and the porewater ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ and O₂ concentrations, respectively. Denitrification is not considered in this model version, but O₂ is not consumed below a threshold, somewhat reflecting the process of denitrification without modeling ${\mathrm{NO}}_{\mathrm{3}}^{-}$. The respective reaction rate parameters for CaCO₃ dissolution and POC oxidation are global constants (see Roth et al., 2014; 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⁻¹ in the form of CaCO₃, 0.24 GtC yr⁻¹ in the form of POC, and 6.96 Tmol Si yr⁻¹ in the form of opal, all within the observational range. The global inventories in the interactive sediment layers amount to 916 GtC for CaCO₃, 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¹³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₃, of CaSiO₃, and to volcanic CO₂ 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 $\mathrm{1}:\mathrm{117}:-\mathrm{17}$ for organic material. Similarly, the silicon (Si) flux is assigned to CaSiO₃ weathering, and the related ALK flux is computed using Si : ALK is equal to 1:2 based on the simplified equation for CaSiO₃ weathering: 2CO₂ + H₂O + CaSiO₃ → Ca²⁺ + $\mathrm{2}{\mathrm{HCO}}_{\mathrm{3}}^{-}$  + SiO₂ (Colbourn et al., 2013). The remaining ALK flux is attributed to CaCO₃ weathering with the stoichiometric ratio C : ALK equal to 1:2 following from CO₂ + H₂O + CaCO₃ → Ca²⁺ + $\mathrm{2}{\mathrm{HCO}}_{\mathrm{3}}^{-}$. 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⁻¹ for weathering of organic material, 0.14 GtC yr⁻¹ for CaCO₃ weathering, 0.08 GtC yr⁻¹ for volcanic CO₂ outgassing, and 6.96 Tmol Si yr⁻¹ for CaSiO₃ weathering. The isotopic signature of the weathering carbon corresponds to the respective signature of the burial fluxes and amounts to δ¹³C = −9.2 ‰, intermediate between isotopically light organic carbon (δ¹³C = −20.5 ‰) and heavier CaCO₃ (δ¹³C = 2.9 ‰). During experiments, weathering fluxes of CaCO₃ and CaSiO₃, and accordingly DIC, DI¹³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, 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 their bioclimatic limits for resources. CO₂ 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 pattern-scaling 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).

2.2 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₂ is prescribed to 284.7 ppm, δ¹³C to −6.305 ‰, and Δ¹⁴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₂ 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 Overview of simulations.

^a The number of ¹⁴C atoms corresponding to 100 GtC with Δ¹⁴C = 0 ‰ is added. ^b Setup without CaCO₃ and CaSiO₃ weathering feedbacks. ^c “Closed system”: atmosphere–ocean–land without ocean sediment module and CaCO₃ and CaSiO₃ weathering feedbacks. ^d Setup with four-box land biosphere instead of LPX. ^e Setup with four-box land biosphere instead of LPX under LGM conditions.

Table 1 summarizes run names and key characteristics of all simulations. Pulse sizes are varied between −250 GtC (removal) to +500 GtC (addition) with δ¹³C of −24 ‰ and Δ¹⁴C of 0 ‰. This corresponds to a hypothetical, sudden uptake or release of carbon from the land biosphere with a typical C₃-like δ¹³C signature. For simplicity, and to ease interpretation, we set the Δ¹⁴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 ¹⁴C of “young” plant, litter, and soil material. In sensitivity experiments, Δ¹⁴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 isotopic signatures are set to the signature for the modern mix of fossil fuels. Finally, we also performed a run where only ¹⁴C, but no ¹²C or ¹³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₃ and CaSiO₃ 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” atmosphere–ocean–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₅₀₀ and 4box₅₀₀. In 4box₅₀₀ the Bern3D setup and spin-up is as in WEA₅₀₀, except that an atmospheric CO₂ 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₂ fertilization, instead of LPX-Bern is coupled to the Bern3D in these runs. In LGM₅₀₀ the model is forced by glacial boundary conditions instead of preindustrial conditions as in 4box₅₀₀ and in WEA₅₀₀. CO₂ 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₅₀₀ should be treated with some caution. Differences in results between WEA₅₀₀ and 4box₅₀₀ are due to differences in the two land biosphere models. Differences between 4box₅₀₀ and LGM₅₀₀ are exclusively due to differences in climatic boundary conditions.

2.3 Data reduction

The remaining fraction of a pulse-like perturbation or IRF at time t after the pulse release is defined for atmospheric CO₂ as

$$\begin{array}{}\text{(1)}& \mathrm{IRF}\left({\mathrm{CO}}_{\mathrm{2},\mathrm{a}}\right)\left(t\right)={\displaystyle \frac{\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2},\mathrm{a}}\left(t\right)}{{\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2},\mathrm{a}}}^{\text{ini}}}}={\displaystyle \frac{{\mathrm{CO}}_{\mathrm{2},\mathrm{a}}\left(t\right)-{\mathrm{CO}}_{\mathrm{2},\mathrm{a}}^{\mathrm{ctrl}}\left(t\right)}{P/\left(\mathrm{2.12}\mathrm{GtC}{\mathrm{ppm}}^{-\mathrm{1}}\right)}}.\end{array}$$

Δ 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₂ perturbation in parts per million (ppm) when divided by 2.12 GtC ppm⁻¹. 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 δ¹³C_a is

$$\begin{array}{}\text{(2)}& \mathrm{IRF}\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a}}\right)\left(t\right)={\displaystyle \frac{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a}}\left(t\right)-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a}}^{\text{ctrl}}\left(t\right)}{\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a}}^{\text{ini}}}}.\end{array}$$

$\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a}}^{\text{ini}}$ denotes the initial isotopic perturbation in the atmosphere and equals:

$$\begin{array}{}\text{(3)}& {\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{a}}^{\text{ini}}={\displaystyle \frac{P}{N_{\mathrm{a},\mathrm{0}}+P}}\cdot ({{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_P-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{a},\mathrm{0}}).\end{array}$$

δ¹³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 Δ¹⁴C.

An analytical expression for the IRF is calculated from experiment WEA₅₀₀. The IRF is fitted by the sum of five exponential terms to

$$\begin{array}{}\text{(4)}& \mathrm{IRF}\left({\mathrm{CO}}_{\mathrm{2},a}\right(t\left)\right)=a_{\mathrm{0}}+\sum _{i=\mathrm{1}}^{\mathrm{5}}{a}_i\cdot \exp \left({\displaystyle \frac{-t}{{\mathit{\tau }}_i}}\right)\text{for}t\ge \mathrm{0}\end{array}$$

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₀ is the fraction of the perturbation remaining constant in the atmosphere (τ₀=∞). 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 δ¹³C_a and Δ¹⁴C_a is fitted with n=4 exponents.

Principal component empirical orthogonal function (PC-EOF) analysis is used to fit the ¹³C isotopic perturbations of dissolved inorganic carbon (DIC) in the ocean (${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{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 ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{DIC}}$(t, x). PCs and EOFs are then used to build a cost-efficient substitute model of the spatiotemporal evolution of ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{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.

3 Results

3.1 Earth system response to a pulse release of carbon

Impulse response functions for atmospheric CO₂ and δ¹³C

We first describe the IRF for the atmospheric CO₂ perturbation using experiment WEA₂₅₀ (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.

Figure 1 (a) CO₂, (b) δ¹³C, and (c) Δ¹⁴C perturbations in the atmosphere for different simulations with different pulse sizes and model configurations. The isotopic signature is δ¹³C = −24 ‰ for all pulse sizes except the 5000 GtC pulse with δ¹³C = −28 ‰. See 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₅₀₀, SED₅₀₀, CLO₅₀₀, 4box₅₀₀, and LGM₅₀₀) for the time interval 1 to 100 kyr are shown in Fig. 2 for better visibility.

CaCO₃ 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₂ perturbation remain airborne (Fig. 1a). The fit by exponentials shows that a fraction of 11.2 % of the initial perturbation is removed with a timescale of 5.5 kyr. We associate these values with the process of CaCO₃ 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₂ is less than 1 ppm and less than 0.03 ‰ in δ¹³CO₂. 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₅₀₀₀) than results presented by Colbourn et al. (2015) and Lord et al. (2016). The remaining airborne fraction of the CO₂ perturbation is 5 % and less after 100 kyr in all studies.

Compared to CO₂, the atmospheric δ¹³C perturbation is initially removed much faster (Fig. 1b). Within the first decades, IRF(δ¹³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(${{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{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.

Table 2 Fit parameters for IRF(ΔCO₂(t)) and IRF(Δδ¹³C_a(t)) of WEA₅₀₀, IRF(ΔΔ¹⁴C_a,dead(t)), and IRF(ΔΔ¹⁴C_a,only(t)). Note that in the case of IRF(ΔΔ¹⁴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).

Processes

Next, we address the processes behind the different temporal evolution of CO_2,a and δ¹³C_a. We focus on the ocean as ocean uptake is thought to dominate on long timescales. A perturbation in atmospheric CO₂ or in ¹³CO₂ 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 δ¹³C_a is linked to the aquatic carbonate chemistry and the associated equilibrium between dissolved CO₂ and bicarbonate and carbonate ions (CO₂ + H₂O ⇔ ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ + H⁺ ⇔ ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ + 2H⁺) (Dickson et al., 2007). The removal of an atmospheric CO₂ 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 ¹³CO₂∕¹²CO₂ is hardly affected by this chemical buffering, because ¹³CO₂ and ¹²CO₂ are affected approximately equally by the acid–base chemical reactions. The chemical buffering diminishes the uptake capacity of the ocean for CO₂ but not for the isotopic ratio. Correspondingly, the IRF for CO_2,a decreases less rapidly than the IRF for δ¹³C_a, and IRF(δ¹³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 δ¹³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₂. For CO₂, the expression reads as

$$\begin{array}{}\text{(5)}& {\text{IRF}}_{\mathrm{\infty }}\left({\mathrm{CO}}_{\mathrm{2},\mathrm{a}}\right)={\displaystyle \frac{N_{\mathrm{a},\mathrm{0}}}{N_{\mathrm{a},\mathrm{0}}+\frac{\mathrm{1}}{\mathit{\xi }}{N}_{\mathrm{o},\mathrm{0}}}}\text{for}\mathrm{0}\le t\le \mathrm{2}\text{kyr}.\end{array}$$

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₂ partial pressure and in dissolved inorganic carbon (DIC) (Revelle and Suess, 1957). ξ varies with environmental conditions and increases with the size of the CO₂ perturbation. ξ is on the order of 10 for small pulse sizes (here −250 to 500 GtC). The dependency of IRF(CO₂) 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_{\mathrm{o},\mathrm{0}}/\mathit{\xi }$). The smaller this ratio, the larger IRF(CO_2,a) is.

For δ¹³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):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& {\text{IRF}}_{\mathrm{\infty }}\left({{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{a}}\right)={\displaystyle \frac{N_{\mathrm{a},\mathrm{0}}+P}{N_{\mathrm{a},\mathrm{0}}+N_{\mathrm{o},\mathrm{0}}+N_{\mathrm{b},\mathrm{0}}+P}}\\ & \cdot {\displaystyle \frac{({{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_P-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\text{mean}})}{({{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_P-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{a},\mathrm{0}})}},\text{for}\mathrm{0}\le t\le \mathrm{2}\text{kyr}.\end{array}\end{array}$$

Equation (6) shows that IRF(δ¹³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. δ¹³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₂ and ¹³CO₂ are both about equally affected by the aquatic chemistry. In turn, their ratio, i.e., δ¹³C, is hardly affected by this acid–base buffering.

Equations (5) and (6) reveal the fundamental difference between a perturbation in CO₂ (and ¹³CO₂) and in the isotopic ratio (δ¹³CO₂). For small pulse sizes (P→0), neglecting the influence of the land biosphere and the correction term in Eq. (6) for δ¹³C, Eqs. (6) and (5) are formally equal, except that a buffer factor ξ of 1, instead of ∼12, applies for δ¹³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 in the case of δ¹³C_a but proportionally to N_a,0 and $N_{\mathrm{o},\mathrm{0}}/\mathit{\xi }$ in the case of CO₂. As a consequence, the perturbation in the ratio δ¹³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₂ 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₅₀₀ versus WEA₅₀₀₀ in Fig. 1a after 1–2 kyr). Regarding δ¹³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 $\le t\le$ 20 kyr), CaCO₃ compensation further reduces the atmospheric CO₂ perturbation and removes the differences in IRF(CO_2,a) arising from different pulse sizes (Fig. 1a). The process of CaCO₃ compensation is briefly explained as follows. CO₂ is taken up by the ocean and partly reacts to form bicarbonate (${\mathrm{HCO}}_{\mathrm{3}}^{-}$) and hydrogen (H⁺) ions. This makes the water more acidic; the carbonate ion concentration ([${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$]) is reduced and the water becomes corrosive to CaCO₃. In turn, sedimentary CaCO₃ dissolves to ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ and Ca²⁺ ions, partly removing the perturbations in [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] and [H⁺], a so-called seafloor neutralization. Under more acidic conditions, less alkalinity and carbon, in the form of CaCO₃, 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 [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] and [H⁺] , a so-called 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₂ 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₃ compensation for the pulse perturbation is completed (see Appendix):

$$\begin{array}{}\text{(7)}& {\text{IRF}}_{{\mathrm{CaCO}}_{\mathrm{3}}}\left({\mathrm{CO}}_{\mathrm{2},\mathrm{a}}\right)=\mathrm{2}\cdot \mathrm{1.5}\cdot {\displaystyle \frac{N_{\mathrm{a},\mathrm{0}}}{N_{\mathrm{o},\mathrm{0}}}},\text{for}t\sim \mathrm{20}\text{kyr}.\end{array}$$

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₃, equivalent to $\sim \mathrm{0.5}\times 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₃ compensation depends on the ratio of the initial atmosphere-to-ocean carbon inventory (N_a,0 to N_o,0). Therefore, IRF${}_{{\mathrm{CaCO}}_{\mathrm{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(δ¹³C_a) stays nearly constant between 1 and 10 kyr (Fig. 1b). This implies that CaCO₃ compensation has little influence on the removal of the δ¹³C perturbation. Reasons are that the change in ocean carbon inventory due to CaCO₃ dissolution (burial) is small compared to the total ocean inventory, and that the isotopic signature of CaCO₃ is similar to that of DIC.

On even longer timescales, the CaSiO₃ 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 δ¹³C_a is removed by the burial flux of organic carbon and CaCO₃. 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 δ¹³C (Δδ¹³C) within the ocean as further discussed below in the context of the factorial simulations and in the Appendix.

Figure 2 (a) CO₂ and (b) δ¹³C perturbations in the atmosphere for the factorial simulations WEA₅₀₀, SED₅₀₀, CLO₅₀₀, 4box₅₀₀, and LGM₅₀₀. 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₅₀₀ was run for 90 kyr instead of 100 kyr.

Factorial simulations

We now turn to the factorial simulations to further quantify the effect of CaCO₃ compensation, burial fluxes, and CaCO₃ and CaSiO₃ weathering feedbacks on the CO₂ and δ¹³C perturbations. The evolution of IRF(CO_2,a) for experiment SED₅₀₀ (no terrestrial weathering feedback but with marine sediments) is comparable to WEA₅₀₀ until simulation year 1000 and, it differs thereafter due to varying weathering rates in WEA₅₀₀ (Figs. 1a and 2a). There is no addition of alkalinity to the ocean from enhanced terrestrial weathering in SED₅₀₀, and ΔCO_2,a levels out at a higher value than in WEA₅₀₀, with about 8 % of the perturbation remaining constantly airborne in SED₅₀₀. 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₃ compensation leads, in both WEA₅₀₀ and SED₅₀₀, 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₃ 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₅₀₀ 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 ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{atm}}$ of SED₅₀₀ and WEA₅₀₀, the evolution is initially comparable, but on multimillennial timescales the isotopic perturbation is, somewhat surprisingly and in contrast to the CO₂ perturbation, removed slower in WEA₅₀₀ compared to SED₅₀₀ (Fig. 2b). In the case of SED₅₀₀ 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⁻¹). In simulation WEA₅₀₀, with weathering feedbacks enabled, the long-term removal timescale of ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{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₅₀₀. In the case of SED₅₀₀ 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 δ¹³C_a for these small pulse sizes and with δ¹³C = −24 ‰ is rather small after a few thousand years. For example, in the case of SED₅₀₀, ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{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₅₀₀, 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₂ perturbation for experiment CLO₅₀₀ 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). δ¹³C_a levels out around 2 kyr at around 3 %, which is as expected from Eq. (6). The further decrease on multimillennial timescales results from model drift in the LPX-Bern model.

Next we compare WEA₅₀₀, where LPX-Bern is used as land model, with 4box₅₀₀, where a four-box land biosphere is used. Differences in IRF(CO_2,a) between WEA₅₀₀ and 4box₅₀₀ are due to the different land biosphere models. The four-box land biosphere model takes up less of the CO₂ 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(δ¹³C_a) (Fig. 2b) are largest in the first decades after the pulse, and the root-mean-square deviation in IRF(δ¹³C_a) amounts to less than 1 % between the two simulations.

Finally, we compare the simulations LGM₅₀₀, forced by glacial boundary conditions, and 4box₅₀₀, forced by preindustrial conditions. Differences between the two simulations are solely due to the difference in climate and CO₂ forcing. IRF(CO_2,a) is lower by a few percent in LGM₅₀₀ compared to 4box₅₀₀ (Fig. 2a). The difference is largest in the first 10 kyr, and it is up to 12 % of the perturbation. The root-mean-square deviation between the two simulations amounts to ∼2 %. The larger oceanic carbon uptake in LGM₅₀₀ compared to 4box₅₀₀ 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(δ¹³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_∞(δ¹³C_a) = 0.027 compared to IRF_∞(δ¹³C_a) = 0.032 when using the preindustrial value (589.4 GtC) in the equation. These estimates of IRF(δ¹³C_a) are in good agreement with the results from LGM₅₀₀ and 4box₅₀₀ (Fig. 2b). In summary, the influence of glacial compared to preindustrial boundary conditions on IRF(δ¹³C_a) appears modest in our model.

Figure 3 Inventories of the (a) carbon and (b) δ¹³C perturbations for experiment WEA₂₅₀. Shown are inventories of the perturbations in the atmosphere (light blue), the land biosphere (green), and the ocean (blue). Hatching indicates the influence of changes in the lithosphere (including sediments) due to CaCO₃ compensation and imbalances in the weathering and burial fluxes on the carbon and isotopic inventories of the ocean–atmosphere–land system. For carbon, this inventory is elevated above the pulse perturbation of 250 GtC until about 70 kyr and depleted thereafter. For δ¹³C, the atmosphere–ocean–land inventory is continuously reduced over the course of the simulation. Magenta lines show results for experiment WEA₋₂₅₀ multiplied by −1.

The budgets of the carbon and carbon isotope perturbations

Next, we discuss the budgets of the carbon and ¹³C perturbations. Figure 3 shows results for simulation WEA₂₅₀, with a pulse release of 250 GtC, and for simulation WEA₋₂₅₀, 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₂₅₀. 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 compensation and enhanced terrestrial carbonate and silicate weathering increase ocean alkalinity, and thereby the uptake capacity for atmospheric CO₂ and the remaining atmospheric CO₂ 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₃ 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 inventory 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 δ¹³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–ocean–land system. The δ¹³C perturbation of the burial flux (combined POC and CaCO₃) follows the negative δ¹³C perturbation in surface DIC. Therefore, the burial flux is what ultimately removes the δ¹³C perturbation from the atmosphere–ocean–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₅₀₀) and without weathering feedbacks (SED₅₀₀) in Fig. 4. The inventory in the relatively fast-exchanging 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₅₀₀, 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₃. In contrast, CaCO₃ burial is enhanced in simulation WEA₅₀₀, 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₅₀₀ (∼10 000 GtC ‰) and WEA₅₀₀ (∼7000 GtC ‰). Changes in CaCO₃ 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 fast-exchanging reservoirs at 100 kyr in simulations SED₅₀₀ and WEA₅₀₀.

The contributions of the anomalies in the POC and CaCO₃ burial fluxes may be approximately attributed to changes in the signature (Δδ) and to changes in flux (ΔF) relative to the initial flux, F₀, and signature, δ₀, before the perturbation:

$$\begin{array}{}\text{(8)}& \mathrm{\Delta }(F\cdot \mathit{\delta })=F_{\mathrm{0}}\cdot \mathrm{\Delta }\mathit{\delta }+\mathrm{\Delta }F\cdot {\mathit{\delta }}_{\mathrm{0}}.\end{array}$$

The isotopic perturbation of the POC burial flux is mainly mediated by a change in the signature of the organic carbon buried (Δδ¹³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₅₀₀. 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₃ burial flux and changes in CaCO₃ burial contribute significantly to the isotopic perturbation. At the end of simulation WEA₅₀₀, the mean δ¹³C signature of CaCO₃ burial is on average 0.14 ‰ more negative than the input from weathering. In turn, the term F₀⋅Δδ for CaCO₃ burial contributes $\sim -\mathrm{3100}$ GtC ‰ to the mitigation of the initial isotopic perturbation. This is partly counteracted by excess burial of CaCO₃, with ΔF⋅δ₀ equal to ∼1300 GtC ‰. This yields an overall burial contribution of $\sim -\mathrm{1800}$ GtC ‰. The δ¹³C perturbation is further removed by increased weathering input of CaCO₃ (see “Model and experimental setup”), contributing about 1100 GtC ‰ over 100 years in WEA₅₀₀.

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₃ cycle to the isotopic budget are related to both a change in the signature of the mean CaCO₃ burial flux and to changes in burial and weathering carbon fluxes.

Figure 4 The perturbation budget for (a) carbon and (b) carbon isotopes. Solid lines show results from experiment WEA₅₀₀; dashed lines are from experiment SED₅₀₀. (a, b) The pulse (500, −12 000 GtC ‰) is subtracted from the combined atmosphere–ocean–land-reactive sediment inventory. Burial fluxes are plotted inversely, and the sum of burial (green, orange lines) and weathering (black) fluxes yields the combined inventory (pink).

Influence of pulse size

The IRF for atmospheric CO₂ and δ¹³C is sensitive to the magnitude of the pulse size (Fig. 1a, b), and we next address related differences between simulations. For CO₂, 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₅₀₀₀) 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₅₀₀₀ compared to smaller pulse sizes (Fig. 1a).

Second, productivity on land depends nonlinearly on CO₂ (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₂₅₀ and WEA₋₂₅₀ 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₂ concentrations, and thus air–sea gas exchange, and partly leveling out the effect from the land biosphere.

For δ¹³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(δ¹³C_a) of ∼2 % for the pulse addition of 250 GtC (WEA₂₅₀) but of only ∼1 % for a corresponding removal (WEA₋₂₅₀).

Radiocarbon

Finally, we discuss the Δ¹⁴C perturbation experiments. In experiment ¹⁴C_only, a positive pulse of ¹⁴C is added to the atmosphere, but there is no perturbation in CO₂, analogous to a radiocarbon production pulse by atomic bomb tests or cosmic rays. The initial evolution of the IRF(Δ¹⁴C_a) is very similar to the one of IRF(δ¹³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 ¹⁴C₋₅₀₀ and ¹⁴C_dead, a radiocarbon-depleted carbon source of 100 GtC is added to the atmosphere causing a negative perturbation in Δ¹⁴C and a positive perturbation in CO₂ in the atmosphere. As visible from Fig. 1a for WEA₁₀₀, a small percentage of the CO₂ perturbation remains airborne even after 100 kyr. This explains the persistence of a small atmospheric perturbation in Δ¹⁴C, well beyond the lifetime of the initially added ¹⁴C. IRF(Δ¹⁴C_a) is larger for ¹⁴C₋₅₀₀ than for ¹⁴C_dead on these very long timescales. This is because the perturbation in CO₂, and in turn the long-term perturbation in Δ¹⁴C_a, is the same for both experiments. However, the initial perturbation in Δ¹⁴C – the denominator of the response function (see Eq. 2) – is smaller in ¹⁴C₋₅₀₀ than in ¹⁴C_dead.

3.2 The δ¹³C perturbation of marine dissolved inorganic carbon and its approximation by PC-EOF

The ocean mean perturbation in δ¹³C of DIC (${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{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₅₀₀) occurs ∼1500 years after the pulse release (Fig. 5), reflecting ocean mixing timescales. Thereafter, the δ¹³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 δ¹³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 δ¹³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, ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{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 ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{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 air–sea gas exchange at higher temperatures, leading to more negative δ¹³C signatures in the surface ocean. Additionally, gross air–sea fluxes are increased as a result of the higher CO₂ concentrations, decreasing the air–sea disequilibrium. This leads to more negative δ¹³C signatures in mid- and low-latitude surface waters and less negative signatures in cold high-latitude waters (see Menviel et al., 2015). This surface-to-deep ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{DIC}}$ gradient is dampened in the case of experiment WEA₅₀₀ by enhanced weathering of CaCO₃ on land, adding carbon with a positive δ¹³C signature (∼2.9 ‰) to the surface ocean. As discussed above, these spatial Δδ¹³C gradients are important and co-govern the removal rate of the δ¹³C perturbation. The larger-than-average δ¹³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₂ concentrations is higher, leading to lower δ¹³C values of the POC export flux. This accelerates the removal of the δ¹³C perturbation in comparison to a uniformly mixed ocean.

Figure 5 Modeled global-average perturbation of the isotopic signature of dissolved inorganic carbon in the ocean (δ¹³C_DIC; thick gray line) in response to a carbon pulse to the atmosphere of 500 GtC with a δ¹³C signature of −24 ‰. The thin colored lines show the perturbation as reconstructed with one (blue), two (green), and three (orange) principal components.

Next, the temporal evolution of the δ¹³C_DIC perturbation and its spatial pattern are described by PC-EOF (Fig. 7). The first PC-EOF pair captures approximately the δ¹³C_DIC evolution after 1 kyr, while the second and third PC-EOF pairs capture, together with the first pair, the decadal-to-centennial-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 ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{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 δ¹³C values in surface and intermediate waters in the first centuries of the simulation, whereas on longer timescales the depth gradient in ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{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).

Figure 6 Modeled δ¹³C_DIC perturbations (‰) (a, c, e, g) and difference between modeled and reconstructed δ¹³C_DIC perturbation using the first three principal components. These differences are shown in percent of the modeled perturbation (b, d, f, h). Shown are time slices at 100 years (a, b), 1 kyr (c, d), 10 kyr (e, f), and 50 kyr (g, h) after the pulse release of carbon. Note that in panel (b) perturbation values smaller than $\left|\mathrm{0.05}\right|$ ‰ (abyssal Atlantic and Pacific) were masked for better visibility as the perturbation has not yet propagated into these regions (see panel a).

The global-mean ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{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-mean-square error (RMSE) between reconstructed and modeled global-mean ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{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 ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{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 ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{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 ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{DIC}}$ signal in these water masses. The performance of the reconstruction increases over the course of the simulation. The RMSE of the ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{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 ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{DIC}}$ with three PC-EOF pairs shows good results, and deviations from the modeled ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{DIC}}$ are of the order of 10 % of the perturbation and smaller on centennial-to-millennial timescales (Fig. 6d, f, h). To reconstruct the modeled evolution of ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{DIC}}$ precisely during the first few centuries after a perturbation, more than three PC-EOF pairs are necessary. In summary, the spatiotemporal evolution of ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{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.

Figure 7 EOF pattern of the (a) first, (c) second, and (d) third principal components used to represent the spatial perturbation in Δδ¹³C of dissolved inorganic carbon in the ocean. Corresponding PC time series are shown in (b).

4 Discussion and conclusion

The aim of this study is to investigate the evolution of a δ¹³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 δ¹³C response is compared to that of atmospheric CO₂ and radiocarbon.

The impulse response (or Green's) functions (IRF) for atmospheric CO₂ and its δ¹³C signature are fitted by a sum of exponential functions. Additionally, we show that the spatiotemporal evolution of a δ¹³C_DIC perturbation is reasonably represented with three PC-EOF pairs. Deviations between the PC-EOF reconstruction and modeled δ¹³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 δ¹³C_DIC variations from marine sediments in future studies. The evolution of the atmospheric CO₂ 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., 1998, 2009; Ridgwell 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 (Joos et al., 1996).

We quantify the processes leading to the different removal timescales of atmospheric CO₂ and δ¹³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 ¹³C∕¹²C is not buffered by marine carbonate chemistry, unlike CO₂. The atmospheric δ¹³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₂ perturbation remain airborne for millennia. On timescales of a few millennia to 20 kyr, the initial atmospheric CO₂ perturbation is lowered by carbonate compensation to a fraction of about 8 %. In contrast, this chemical buffering hardly affects the atmospheric δ¹³C perturbation. On even longer timescales, the CO₂ and the isotopic perturbation are removed by different, though related, processes. For δ¹³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₂, the excess weathering of silicate rocks in response to the CO₂ and climate perturbation adds alkalinity to the ocean, which leads to a complete removal of the atmospheric CO₂ perturbation.

Gradients in the δ¹³C_DIC perturbation strongly influence the long-term removal timescale of the isotopic perturbation. In the pulse experiments, Δδ¹³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₂. Changes in CaCO₃ weathering, on the other hand, diminish the surface perturbation (e.g., experiment WEA₅₀₀ vs. SED₅₀₀ in Fig. 1a). All these processes lead to spatial gradients in the δ¹³C perturbation in the ocean. As a result of these gradients in Δδ¹³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 δ¹³C signature of the relentless flux of organic and calcium carbonate particles will influence the δ¹³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 δ¹³C between similar climate states such as the Penultimate Glacial Maximum (PGM) and the Last Glacial Maximum (LGM). Substantial temporal δ¹³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 δ¹³C signature of weathering and burial fluxes, varying contribution of volcanic outgassing of CO₂, and changes in the amount of carbon stored in the land biosphere, especially in yedoma and permafrost soils – have been discussed for the PGM–LGM δ¹³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 δ¹³C would require an opposing δ¹³C change in the ocean. This appears to be in conflict with marine and ice core records, which suggest that δ¹³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 δ¹³C in the surface ocean, and in turn δ¹³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 δ¹³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 δ¹³C difference. Schneider et al. (2013) estimated required changes in land carbon storage to match the ice core and marine temporal offsets in δ¹³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 δ¹³C offset (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 δ¹³C differences between PGM and LGM. However, further work, which considers the transient evolution of CO₂, δ¹³C, and other proxies, is required to gain further insight into the contributions of individual mechanisms to long-term δ¹³C changes recorded in ice and marine cores.

Further, the different behavior of CO₂ and δ¹³C_atm perturbations also has consequences for centennial-scale CO₂ and δ¹³C_atm variations such as during Heinrich stadials 4 and 1 (HS4 and HS1). Variations in CO₂ and δ¹³C_atm during HS4 and HS1 have been measured on Antarctic ice cores. The data show an increase in CO₂ on the order of ∼10 ppm over 200–300 years, accompanied by a decrease in δ¹³C_atm of ∼0.2 ‰. This results in a value of roughly −0.02 ‰ ppm⁻¹ for the ratio of the change in δ¹³C to the change in CO₂ ($r=\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}/\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}$). A release of terrestrial carbon has been discussed as a possible cause of these events (Bauska et al., 2016, 2018).

Figure 8 (a) Temporal evolution of the ratio of Δδ¹³C to ΔCO₂ in the atmosphere. Results are from simulations WEA₁₀₀ forced by preindustrial boundary conditions (solid lines) and LGM₁₀₀ forced by Last Glacial Maximum boundary conditions (thin dashed lines). 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 Δδ¹³C-to-ΔCO₂ ratio of −0.2 ‰ to 10 ppm as seen in ice core records for HS1 and HS4. (b) Response in ΔCO₂ and ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{atm}}$ to a uniform release of 20, 40, 80, 120, and 160 GtC (colors) of isotopically light carbon (δ¹³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. Δδ¹³C-to-ΔCO₂ change as seen in ice core records for HS1 and HS4 (Bauska et al., 2016, 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 pulse-like 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⁻¹ 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 ¹³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⁻¹ thereafter. Thus, r recorded in ice after a hypothetical pulse release of carbon would be above −0.02 ‰ ppm⁻¹ 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={\int }_{\mathrm{0}}^T\mathrm{d}t\cdot e\left(t\right)\cdot$ IRF(δ¹³C)$(T-t)/{\int }_{\mathrm{0}}^T\mathrm{d}t\cdot e\left(t\right)\cdot$ IRF(CO₂)(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⁻¹ 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 terrestrial carbon with δ¹³C $=-\mathrm{24}$ ‰ over 100 to 400 years yield a linear response in ΔCO₂ and ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{atm}}$ (Fig. 8b) with a slope r of about −0.01 ‰ ppm⁻¹, 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 δ¹³C depletion in vegetation and soils may vary between −24 ‰ for heavily discriminated C₃ plant material and to −5 ‰, which is indicative of C₄ plant material (Keller et al., 2017). C₄ 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 δ¹³C_atm more than in CO₂ (see Köhler et al., 2011). On the other hand, the CO₂ 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₂ events or not. Other mechanisms in addition to a terrestrial release may have contributed to the CO₂ and δ¹³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₂ and δ¹³C are similar under preindustrial and glacial CO₂ in our model. The direct influence of lower atmospheric CO₂ on IRF(δ¹³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₂ and of δ¹³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 δ¹³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 δ¹³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 δ¹³C during glacial–interglacial, and longer, periods.

Appendix A: Analytical expressions for the impulse response functions of CO_2,a and δ¹³C_a

Ocean invasion of an atmospheric CO₂ perturbation

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₂ (Eq. A7) and δ¹³C (Eq. A29), starting in this section with IRF(CO₂).

For small pulses, changes in sea surface temperature, in CO₂ 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.

Figure A1 Hovmöller-type diagram showing the temporal evolution of the global-mean vertical profile of ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{DIC}}$ for experiments (a) WEA₅₀₀ and (b) SED₅₀₀.

We consider a two-box model of the atmosphere (index: a) and ocean (index: o) with area $A=A_{\mathrm{a}}=A_{\mathrm{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₂ partial pressure is $p{\mathrm{CO}}_{\mathrm{2},\mathrm{0}}=p{\mathrm{CO}}_{\mathrm{2},\mathrm{a},\mathrm{0}}=p{\mathrm{CO}}_{\mathrm{2},\mathrm{o},\mathrm{0}}$. The initial concentration of dissolved inorganic carbon (DIC) in the ocean is DIC₀, and the carbon inventory in the ocean is given by $N_{\mathrm{o}}=\mathrm{DIC}\cdot h_{\mathrm{o}}\cdot A$. We relate the atmospheric carbon inventory, N_a, to a time-varying atmospheric concentration C_a and to the fixed volume of the model atmosphere: $N_{\mathrm{a}}=C_{\mathrm{a}}\cdot h_{\mathrm{a}}\cdot A$. We are free to select the units of the model-specific concentration C_a. Here, we express C_a in units of DIC₀ and set the initial atmospheric concentration C_a,0 equal to DIC₀. 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:

$$\begin{array}{}\text{(A1)}& C_{\mathrm{a}}={\displaystyle \frac{p{\mathrm{CO}}_{\mathrm{2},\mathrm{a}}}{p{\mathrm{CO}}_{\mathrm{2},\mathrm{0}}}}{\mathrm{DIC}}_{\mathrm{0}},\end{array}$$

where C_a is specifically defined for our box model and should not be confused with the concentration or mixing ratio of CO₂ in the real-world atmosphere. The atmospheric scale-height h_a is given by the requirement that volume times concentration equals inventory:

$$\begin{array}{}\text{(A2)}& h_{\mathrm{a}}={\displaystyle \frac{N_{\mathrm{a},\mathrm{0}}}{A\cdot {\mathrm{DIC}}_{\mathrm{0}}}}.\end{array}$$

The perturbation (Δ) from the initial equilibrium in pCO₂, and similarly for other quantities, is defined by ΔpCO₂(t) = pCO₂(t) − pCO₂(t₀). The perturbation in pCO_2,o is related to that in DIC by the Revelle factor ξ (Revelle and Suess, 1957):

$$\begin{array}{}\text{(A3)}& \mathit{\xi }={\displaystyle \frac{\mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2},\mathrm{o}}}{p{\mathrm{CO}}_{\mathrm{2},\mathrm{0}}}}\cdot {\displaystyle \frac{{\mathrm{DIC}}_{\mathrm{0}}}{\mathrm{\Delta }\mathrm{DIC}}}.\end{array}$$

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, $\mathrm{\Delta }{f}_{\mathrm{a}\to \mathrm{o},\mathrm{net}}$:

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}& \mathrm{\Delta }{f}_{\mathrm{a}\to \mathrm{o},\mathrm{net}}=k_{\mathrm{g}}\cdot (\mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2},\mathrm{a}}-\mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2},\mathrm{o}})\\ & =g\cdot (\mathrm{\Delta }{\mathrm{C}}_{\mathrm{a}}-\mathit{\xi }\mathrm{\Delta }\mathrm{DIC}),\end{array}\end{array}$$

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:

$$\begin{array}{}\text{(A5)}& g=k_{\mathrm{g}}\cdot {\displaystyle \frac{p{\mathrm{CO}}_{\mathrm{2},\mathrm{0}}}{{\mathrm{DIC}}_{\mathrm{0}}}}.\end{array}$$

The IRF, or equivalent airborne fraction, is defined by

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}& \text{IRF}\left({\mathrm{CO}}_{\mathrm{2},\mathrm{a}}\right(t\left)\right)={\displaystyle \frac{\mathrm{\Delta }{N}_{\mathrm{a}}\left(t\right)}{P}}={\displaystyle \frac{\mathrm{\Delta }{N}_{\mathrm{a}}\left(t\right)}{\mathrm{\Delta }{N}_{\mathrm{a}}\left(t\right)+\mathrm{\Delta }{N}_{\mathrm{o}}\left(t\right)}},\\ & \mathrm{for}t>t_{\mathrm{0}}.\end{array}\end{array}$$

P is the magnitude of the initial carbon pulse at time t₀, 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_{\mathrm{a},\mathrm{0}}={\text{DIC}}_{\mathrm{0}}$ and from Eq. (A4) that $\mathrm{\Delta }{\mathrm{DIC}}_{\mathrm{\infty }}=\mathrm{\Delta }{C}_{\mathrm{a},\mathrm{\infty }}/\mathit{\xi }$. Finally, using this, as well as $\mathrm{\Delta }{N}_{\mathrm{a}}=\mathrm{\Delta }{C}_{\mathrm{a}}\cdot h_{\mathrm{a}}\cdot A$ and $\mathrm{\Delta }{N}_{\mathrm{o}}=\mathrm{\Delta }\mathrm{DIC}\cdot h_{\mathrm{o}}\cdot A=\mathrm{\Delta }{C}_{\mathrm{a}}/\mathit{\xi }\cdot h_{\mathrm{o}}\cdot A$, yields the following for the IRF in an ocean–atmosphere system at equilibrium:

$$\begin{array}{}\text{(A7)}& {\mathrm{IRF}}_{\mathrm{\infty }}\left({\mathrm{CO}}_{\mathrm{2},\mathrm{a}}\right)={\displaystyle \frac{h_{\mathrm{a}}}{h_{\mathrm{a}}+\frac{\mathrm{1}}{\mathit{\xi }}{h}_{\mathrm{o}}}}={\displaystyle \frac{N_{\mathrm{a},\mathrm{0}}}{N_{\mathrm{a},\mathrm{0}}+\frac{\mathrm{1}}{\mathit{\xi }}{N}_{\mathrm{o},\mathrm{0}}}}.\end{array}$$

Distribution of a δ¹³C perturbation within the atmosphere–ocean–land-biosphere system

We now turn to the isotope ¹³C. The isotopic ratio is given by the following (see Mook, 1986):

$$\begin{array}{}\text{(A8)}& {}^{\mathrm{13}}R\left(\mathrm{C}\right){=}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}.\end{array}$$

And the isotopic signature in per mill is given by

$$\begin{array}{}\text{(A9)}& {\mathit{\delta }}^{\mathrm{13}}\mathrm{C}=\left({\displaystyle \frac{{}^{\mathrm{13}}R}{{}^{\mathrm{13}}{R}_{\mathrm{SD}}}}-\mathrm{1}\right)\cdot \mathrm{1000}.\end{array}$$

¹³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

$$\begin{array}{}\text{(A10)}& \mathit{\epsilon }=(\mathit{\alpha }-\mathrm{1})\cdot \mathrm{1000}.\end{array}$$

In the following, we neglect for simplicity the approximately 1 % difference between the concentration of ¹²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 (δ¹³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 ¹³C budgets for the three-box model are described by

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{N}_{\mathrm{o}}^i}{\mathrm{d}t}}=F_{\mathrm{a}\to \mathrm{o}}^i-F_{\mathrm{o}\to \mathrm{a}}^i=F_{\mathrm{a}\to \mathrm{o},\mathrm{net}}^i\\ & {\displaystyle \frac{\mathrm{d}{N}_{\mathrm{b}}^i}{\mathrm{d}t}}=F_{\mathrm{a}\to \mathrm{b}}^i-F_{\mathrm{b}\to \mathrm{a}}^i=F_{\mathrm{a}\to \mathrm{b},\mathrm{net}}^i\\ & {\displaystyle \frac{\mathrm{d}{N}_{\mathrm{a}}^i}{\mathrm{d}t}}=-F_{\mathrm{a}\to \mathrm{o},\mathrm{net}}^i-F_{\mathrm{a}\to \mathrm{b},\mathrm{net}}^i+\mathit{\delta }\left(t\right)\cdot P^i.\end{array}\end{array}$$

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 ¹³C. F denotes fluxes between reservoirs. δ is the Kronecker symbol, and Pⁱ is the pulse released to the atmosphere at time $t=t_{\mathrm{0}}=\mathrm{0}$.

The initial atmospheric δ¹³C perturbation

The IRF(t) for ${{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{a}}$ is given by the perturbation in the isotopic signature, ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{a}}\left(t\right)$, divided by the initial (ini) perturbation, ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{a}}^{\text{ini}}$; ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{a}}^{\text{ini}}$ is the perturbation immediately after a carbon input of amount P and with signature ${}^{\mathrm{13}}{R}_P({{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_P$). Mass balance implies that the amount of ¹³C before and after the pulse release is equal:

$$\begin{array}{}\text{(A12)}& N_{\mathrm{a},\mathrm{0}}\cdot {}^{\mathrm{13}}{R}_{\mathrm{a},\mathrm{0}}+P\cdot {}^{\mathrm{13}}{R}_P=(N_{\mathrm{a},\mathrm{0}}+P)\cdot {}^{\mathrm{13}}{R}_{\mathrm{a}}^{\text{ini}}.\end{array}$$

We convert the ratios to δ units and subtract on both sides $(N_{\mathrm{a},\mathrm{0}}+P)\cdot {\mathit{\delta }}^{\mathrm{13}}$C_a,0 to get the initial perturbation:

$$\begin{array}{}\text{(A13)}& \mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a}}^{\text{ini}}={\displaystyle \frac{P}{N_{\mathrm{a},\mathrm{0}}+P}}\cdot ({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_P-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{0}}).\end{array}$$

${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{a}}^{\text{ini}}$ 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 δ¹³C signature.

The relationship between δ¹³C_a, δ¹³C_o, and δ¹³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.

The gross air-to-sea flux per unit area is

$$\begin{array}{}\text{(A14)}& \begin{array}{rl}& {}^{\mathrm{13}}{f}_{\mathrm{a}\to \mathrm{o}}{=}^{\mathrm{13}}{k}_{\mathrm{g}}{\cdot }^{\mathrm{13}}p{\mathrm{CO}}_{\mathrm{2},\mathrm{a}}=k_{\mathrm{g}}\cdot {\mathrm{CO}}_{\mathrm{2},\mathrm{a}}{\cdot }^{\mathrm{13}}{\mathit{\alpha }}_{\mathrm{a}\to \mathrm{o}}\cdot \\ & {}^{\mathrm{13}}R\left(p{\mathrm{CO}}_{\mathrm{2},\mathrm{a}}\right),\end{array}\end{array}$$

where ${}^{\mathrm{13}}{\mathit{\alpha }}_{\mathrm{a}\to \mathrm{o}}{=}^{\mathrm{13}}{k}_{\mathrm{g}}/k_{\mathrm{g}}$ is the discrimination factor for the gross air-to-sea transfer. The gross sea-to-air flux is

$$\begin{array}{}\text{(A15)}& \begin{array}{rl}& {}^{\mathrm{13}}{f}_{\mathrm{o}\to \mathrm{a}}{=}^{\mathrm{13}}{k}_{\mathrm{g}}{\cdot }^{\mathrm{13}}p{\mathrm{CO}}_{\mathrm{2},\mathrm{o}}\\ & =k_{\mathrm{g}}\cdot p{\mathrm{CO}}_{\mathrm{2},\mathrm{o}}\left({\displaystyle \frac{{}^{\mathrm{13}}{k}_{\mathrm{g}}}{k_{\mathrm{g}}}}{\displaystyle \frac{{}^{\mathrm{13}}R\left(p{\mathrm{CO}}_{\mathrm{2},\mathrm{o}}\right)}{{}^{\mathrm{13}}R\left(\mathrm{DIC}\right)}}\right)\cdot {}^{\mathrm{13}}R\left(\mathrm{DIC}\right)\\ \\ & =k_{\mathrm{g}}\cdot p{\mathrm{CO}}_{\mathrm{2},\mathrm{o}}{\cdot }^{\mathrm{13}}{\mathit{\alpha }}_{\mathrm{o}\to \mathrm{a}}\cdot {}^{\mathrm{13}}R\left(\mathrm{DIC}\right).\end{array}\end{array}$$

The term in parentheses is the discrimination factor for the transfer of carbon from the DIC pool to the atmosphere (¹³α_o→a). At equilibrium (eq), the two gross fluxes cancel each other out. It follows with p${\mathrm{CO}}_{\mathrm{2},\mathrm{a},\mathrm{eq}}$ = p${\mathrm{CO}}_{\mathrm{2},\mathrm{o},\mathrm{eq}}$ and ¹³p${\mathrm{CO}}_{\mathrm{2},\mathrm{a},\mathrm{eq}}$ = ¹³p${\mathrm{CO}}_{\mathrm{2},\mathrm{o},\mathrm{eq}}$ that

$$\begin{array}{}\text{(A16)}& {}^{\mathrm{13}}{R}_{\text{eq}}\left(p{\mathrm{CO}}_{\mathrm{2},\mathrm{a}}\right)={\displaystyle \frac{{}^{\mathrm{13}}{\mathit{\alpha }}_{\mathrm{o}\to \mathrm{a}}}{{}^{\mathrm{13}}{\mathit{\alpha }}_{\mathrm{a}\to \mathrm{o}}}}\cdot {}^{\mathrm{13}}{R}_{\text{eq}}\left(\mathrm{DIC}\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{a},\mathrm{o}}}}\cdot {}^{\mathrm{13}}{R}_{\text{eq}}\left(\mathrm{DIC}\right).\end{array}$$

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

$$\begin{array}{}\text{(A17)}& \begin{array}{rl}& {}^{\mathrm{13}}{F}_{\mathrm{a}\to \mathrm{b}}=F_{\mathrm{a}\to \mathrm{b}}{\cdot }^{\mathrm{13}}{\mathit{\alpha }}_{\mathrm{a}\to \mathrm{b}}\cdot {}^{\mathrm{13}}R\left(p{\mathrm{CO}}_{\mathrm{2},\mathrm{a}}\right)\\ & {}^{\mathrm{13}}{F}_{\mathrm{b}\to \mathrm{a}}=F_{\mathrm{b}\to \mathrm{a}}{\cdot }^{\mathrm{13}}{\mathit{\alpha }}_{\mathrm{b}\to \mathrm{a}}\cdot {}^{\mathrm{13}}R\left(N_{\mathrm{b}}\right).\end{array}\end{array}$$

This yields at equilibrium

$$\begin{array}{}\text{(A18)}& {}^{\mathrm{13}}{R}_{\text{eq}}\left(p{\mathrm{CO}}_{\mathrm{2},\mathrm{a}}\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{a},\mathrm{b}}}}\cdot {}^{\mathrm{13}}{R}_{\text{eq}}\left(N_{\mathrm{b}}\right).\end{array}$$

Equations (A16) and (A18) are readily converted into δ notation:

$$\begin{array}{}\text{(A19)}& \begin{array}{rl}& {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{o},\mathrm{eq}}=\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{a},\mathrm{o}}}{\mathrm{1000}}}\right){\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{eq}}+{\mathit{\epsilon }}_{\mathrm{a},\mathrm{o}}\cong {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{eq}}+{\mathit{\epsilon }}_{\mathrm{a},\mathrm{o}}.\\ & {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{b},\mathrm{eq}}=\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{a},\mathrm{b}}}{\mathrm{1000}}}\right){\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{eq}}+{\mathit{\epsilon }}_{\mathrm{a},\mathrm{b}}\cong {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{eq}}+{\mathit{\epsilon }}_{\mathrm{a},\mathrm{b}}.\end{array}\end{array}$$

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₀) 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:

$$\begin{array}{}\text{(A20)}& \mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{\infty }}=\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{\infty }}\cong \mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{o},\mathrm{\infty }}\cong \mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{b},\mathrm{\infty }}.\end{array}$$

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 Δδ¹³C

The mass balance equations for the ocean–land-biosphere–atmosphere system are given for carbon by

$$\begin{array}{}\text{(A21)}& \begin{array}{rl}& N_{\mathrm{a},\mathrm{0}}+N_{\mathrm{o},\mathrm{0}}+N_{\mathrm{b},\mathrm{0}}+P\\ & =N_{\mathrm{a},\mathrm{\infty }}+N_{\mathrm{o},\mathrm{\infty }}+N_{\mathrm{b},\mathrm{\infty }},\mathrm{for}t>t_{\mathrm{0}};\\ & P=\mathrm{\Delta }{N}_{\mathrm{a}}+\mathrm{\Delta }{N}_{\mathrm{o}}+\mathrm{\Delta }{N}_{\text{b}}.\end{array}\end{array}$$

The second equation describes the perturbation only. We use Eqs. (A16) and (A18) to replace N_o⋅¹³R_o with ${\mathit{\alpha }}_{\mathrm{a},\mathrm{o}}\cdot N_{\mathrm{o}}\cdot {}^{\mathrm{13}}{R}_{\mathrm{a}}$, similarly for the land reservoir, and write the mass balance for ¹³C.

$$\begin{array}{}\text{(A22)}& \begin{array}{rl}& (N_{\mathrm{a},\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{o}}\cdot N_{\mathrm{o},\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{b}}\cdot N_{\mathrm{b},\mathrm{0}})\cdot {}^{\mathrm{13}}{R}_{\mathrm{a},\mathrm{0}}+P\cdot {}^{\mathrm{13}}{R}_P\\ & =(N_{\mathrm{a},\mathrm{\infty }}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{o}}\cdot N_{\mathrm{o},\mathrm{\infty }}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{b}}\cdot N_{\mathrm{b},\mathrm{\infty }})\cdot {}^{\mathrm{13}}{R}_{\mathrm{a},\mathrm{\infty }}\\ & \mathrm{for}t>t_{\mathrm{0}}.\end{array}\end{array}$$

Now we convert to δ units:

$$\begin{array}{}\text{(A23)}& \begin{array}{rl}& (N_{\mathrm{a},\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{o}}\cdot N_{\mathrm{o},\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{b}}\cdot N_{\mathrm{b},\mathrm{0}}+P)\cdot \mathrm{1000}\\ & +(N_{\mathrm{a},\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{o}}\cdot N_{\mathrm{o},\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{b}}\cdot N_{\mathrm{b},\mathrm{0}})\cdot {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{0}}+P\cdot {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_P\\ & =(N_{a,\mathrm{\infty }}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{o}}\cdot N_{\mathrm{o},\mathrm{\infty }}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{b}}\cdot N_{\mathrm{b},\mathrm{\infty }})\cdot \mathrm{1000}\\ & +(N_{\mathrm{a},\mathrm{\infty }}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{o}}\cdot N_{\mathrm{o},\mathrm{\infty }}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{b}}\cdot N_{\mathrm{b},\mathrm{\infty }})\cdot {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{\infty }}.\end{array}\end{array}$$

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.

$$\begin{array}{}\text{(A24)}& \begin{array}{rl}& -\mathrm{\Delta }{N}_{\mathrm{o},\mathrm{\infty }}\cdot \stackrel{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC},\mathrm{0}}-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{0}}}{\overbrace{({\mathit{\alpha }}_{\mathrm{a},\mathrm{o}}-\mathrm{1})\cdot \mathrm{1000}}}\\ & -\mathrm{\Delta }{N}_{\mathrm{b},\mathrm{\infty }}\cdot \stackrel{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{b},\mathrm{0}}-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{0}}}{\overbrace{({\mathit{\alpha }}_{\mathrm{a},\mathrm{b}}-\mathrm{1})\cdot \mathrm{1000}}}\\ & +(N_{\mathrm{a},\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{o}}\cdot N_{\mathrm{o},\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{b}}\cdot N_{\mathrm{b},\mathrm{0}})\\ & \cdot {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{0}}+P\cdot {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_P\\ & =(N_{\mathrm{a},\mathrm{\infty }}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{o}}\cdot N_{\mathrm{o},\mathrm{\infty }}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{b}}\cdot N_{\mathrm{b},\mathrm{\infty }})\cdot {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{\infty }}.\end{array}\end{array}$$

Next, we subtract on both sides the third left-hand side term, $(N_{\mathrm{a},\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{o}}\cdot N_{\mathrm{o},\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{a},\mathrm{b}}\cdot N_{\mathrm{b},\mathrm{0}})\cdot {\mathit{\delta }}^{\mathrm{13}}$C_a,0, as well as the term $P\cdot {\mathit{\delta }}^{\mathrm{13}}{C}_{\mathrm{a},\mathrm{0}}$. We rearrange and use again the mass balance for carbon to get

$$\begin{array}{}\text{(A25)}& \begin{array}{rl}& -\mathrm{\Delta }{N}_{\mathrm{o},\mathrm{\infty }}\cdot ({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC},\mathrm{0}}-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{0}}\cdot (\mathrm{1}+{\mathit{\epsilon }}_{\mathrm{a},\mathrm{o}}/\mathrm{1000}\left)\right)\\ & -\mathrm{\Delta }{N}_{\mathrm{b},\mathrm{\infty }}\cdot ({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{b},\mathrm{0}}-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{0}}\cdot (\mathrm{1}+{\mathit{\epsilon }}_{\mathrm{a},\mathrm{b}}/\mathrm{1000}\left)\right)\\ & -P\cdot ({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a}}-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_P)\\ & =(N_{\mathrm{a},\mathrm{0}}+N_{\mathrm{o},\mathrm{0}}\cdot (\mathrm{1}+{\mathit{\epsilon }}_{\mathrm{a},\mathrm{o}}/\mathrm{1000})+N_{\mathrm{b},\mathrm{0}}\\ & \cdot (\mathrm{1}+{\mathit{\epsilon }}_{\mathrm{a},\mathrm{b}}/\mathrm{1000})+P)\cdot \mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{\infty }}.\end{array}\end{array}$$

The terms with ε∕1000 represent a small correction (<0.02) and can be safely neglected. We replace P with $\mathrm{\Delta }{N}_{\mathrm{a},\mathrm{\infty }}+\mathrm{\Delta }{N}_{\mathrm{o},\mathrm{\infty }}+\mathrm{\Delta }{N}_{\mathrm{b},\mathrm{\infty }}$ in the above equation and exchange the left and right side to get the following mass balance for the δ¹³C perturbation:

$$\begin{array}{}\text{(A26)}& \begin{array}{rl}& (N_{\mathrm{a},\mathrm{0}}+N_{\mathrm{o},\mathrm{0}}+N_{\mathrm{b},\mathrm{0}}+P)\cdot \mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{\infty }}\\ & =\mathrm{\Delta }{N}_{\mathrm{a},\mathrm{\infty }}\cdot \left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_P-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{0}}\right)\\ & +\mathrm{\Delta }{N}_{\mathrm{o},\mathrm{\infty }}\cdot \left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_P-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC},\mathrm{0}}\right)\\ & +\mathrm{\Delta }{N}_{\mathrm{b},\mathrm{\infty }}\cdot \left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_P-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{b},\mathrm{0}}\right)\\ & =P\cdot \left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_P-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\text{mean}}\right),\end{array}\end{array}$$

with

$$\begin{array}{}\text{(A27)}& \begin{array}{rl}& {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\text{mean}}={\displaystyle \frac{\mathrm{\Delta }{N}_{\mathrm{a},\mathrm{\infty }}}{P}}\cdot {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{0}}\\ & +{\displaystyle \frac{\mathrm{\Delta }{N}_{\mathrm{o},\mathrm{\infty }}}{P}}\cdot {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{o},\mathrm{0}}\\ & +{\displaystyle \frac{\mathrm{\Delta }{N}_{\mathrm{b},\mathrm{\infty }}}{P}}\cdot {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{b},\mathrm{0}}.\end{array}\end{array}$$

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 δ¹³C perturbation, Δδ¹³C_∞, is well mixed (Eq., A20) within the total ocean–atmosphere carbon reservoir $(N_{\mathrm{a},\mathrm{0}}+N_{\mathrm{o},\mathrm{0}}+N_{\mathrm{b},\mathrm{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 (δ¹³C_P−δ¹³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 (δ¹³C_P−δ¹³C_mean).

By rearranging Eq. (A26) we get the perturbation in the isotopic signature of the atmosphere–ocean–land system:

$$\begin{array}{}\text{(A28)}& \begin{array}{rl}& \mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{\infty }}={\displaystyle \frac{P}{N_{\mathrm{a},\mathrm{0}}+N_{\mathrm{o},\mathrm{0}}+N_{\mathrm{b},\mathrm{0}}+P}}\\ & \cdot \left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_P-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\text{mean}}\right).\end{array}\end{array}$$

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).

IRF(δ¹³C_a) at equilibrium (t∼2 kyr)

Finally, we use the definition for the IRF and Eqs. (A13) and (A28) to obtain the following expression for the IRF at a new equilibrium for the ocean–atmosphere–land-biosphere system:

$$\begin{array}{}\text{(A29)}& \begin{array}{rl}& {\mathrm{IRF}}_{\mathrm{\infty }}\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a}}\right)\\ & ={\displaystyle \frac{\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{a,\mathrm{\infty }}}{\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a}}^{\text{ini}}}}\\ & ={\displaystyle \frac{N_{\mathrm{a},\mathrm{0}}+P}{N_{\mathrm{a},\mathrm{0}}+N_{\mathrm{o},\mathrm{0}}+N_{\mathrm{b},\mathrm{0}}+P}}\cdot {\displaystyle \frac{({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_P-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\text{mean}})}{({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_P-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a},\mathrm{0}})}}.\end{array}\end{array}$$

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₂ 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 δ¹³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 ¹³C (index: i):

$$\begin{array}{}\text{(A30)}& {\text{IRF}}^i\stackrel{P\to \mathrm{0}}{=}{\displaystyle \frac{N_{\mathrm{a},\mathrm{0}}}{N_{\mathrm{a},\mathrm{0}}+\frac{\mathrm{1}}{{\mathit{\xi }}^i}{N}_{\mathrm{o},\mathrm{0}}}}\cdot M^i.\end{array}$$

However, for δ¹³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 δ¹³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₂ 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⁻³ and an ocean area of 3.6×10¹⁴ m². 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 % of the initial perturbation is still found in the atmosphere after a decade. These numbers are comparable in magnitude to the IRF of CO₂ 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₂ 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 δ¹³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 δ¹³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₂ perturbation

Next, we address CaCO₃ compensation of the carbon added by the pulse. We derive an expression for the IRF(CO_2,a) at the time when CaCO₃ 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₃ compensation, again in equilibrium with the atmosphere. Concentrations of CO₂, ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, and ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ in seawater are related by

$$\begin{array}{}\text{(A31)}& {\displaystyle \frac{p{\mathrm{CO}}_{\mathrm{2},\mathrm{o}}\cdot \left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}{\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]}}={\displaystyle \frac{K_{\mathrm{2}}^{\prime }\cdot K_{\text{H}}}{K_{\mathrm{1}}^{\prime }}}=\text{const}.\end{array}$$

$K_{\mathrm{1}}^{\prime }$ and $K_{\mathrm{2}}^{\prime }$ are the apparent dissociation constants for carbonic acid and K_H is the Henry's law solubility product for CO₂. For simplicity, we assume the ratio of $K_{\mathrm{1}}^{\prime }$, $K_{\mathrm{2}}^{\prime }$, 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

$$\begin{array}{}\text{(A32)}& {\displaystyle \frac{p{\mathrm{CO}}_{\mathrm{2},\mathrm{\infty }}}{p{\mathrm{CO}}_{\mathrm{2},\mathrm{0}}}}\approx {\left({\displaystyle \frac{\left[{\mathrm{HCO}}_{\mathrm{3},\mathrm{\infty }}^{-}\right]}{\left[{\mathrm{HCO}}_{\mathrm{3},\mathrm{0}}^{-}\right]}}\right)}^{\mathrm{2}}\cdot {\displaystyle \frac{\left[{\mathrm{CO}}_{\mathrm{3},\mathrm{0}}^{\mathrm{2}-}\right]}{\left[{\mathrm{CO}}_{\mathrm{3},\mathrm{\infty }}^{\mathrm{2}-}\right]}}.\end{array}$$

Sedimentary dissolution of CaCO₃ in response to the carbon pulse tends to keep the carbonate ion concentration constant over the course of the simulations, yielding

$$\begin{array}{}\text{(A33)}& {\displaystyle \frac{p{\mathrm{CO}}_{\mathrm{2},\mathrm{\infty }}}{p{\mathrm{CO}}_{\mathrm{2},\mathrm{0}}}}\approx {\left({\displaystyle \frac{\left[{\mathrm{HCO}}_{\mathrm{3},\mathrm{\infty }}^{-}\right]}{\left[{\mathrm{HCO}}_{\mathrm{3},\mathrm{0}}^{-}\right]}}\right)}^{\mathrm{2}}.\end{array}$$

${\mathrm{HCO}}_{\mathrm{3}}^{-}$ is by far the dominant form of DIC. This allows one to approximate ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ by DIC in Eq. (A33) to get

$$\begin{array}{}\text{(A34)}& {\displaystyle \frac{p{\mathrm{CO}}_{\mathrm{2},\mathrm{\infty }}}{p{\mathrm{CO}}_{\mathrm{2},\mathrm{0}}}}\approx {\left({\displaystyle \frac{{\mathrm{DIC}}_{\mathrm{o},\mathrm{0}}+\mathrm{\Delta }{\mathrm{DIC}}_{\mathrm{o}}}{{\mathrm{DIC}}_{\mathrm{o},\mathrm{0}}}}\right)}^{\mathrm{2}}={\left({\displaystyle \frac{N_{\mathrm{o},\mathrm{0}}+\mathrm{\Delta }{N}_{\mathrm{o}}}{N_{\mathrm{o},\mathrm{0}}}}\right)}^{\mathrm{2}}.\end{array}$$

Δ refers here to the difference between the final and initial state. With pulse sizes substantially smaller than the oceanic DIC inventory, $\mathrm{\Delta }{N}_{\mathrm{o}}\ll N_{\mathrm{o},\mathrm{0}}$, we approximate

$$\begin{array}{}\text{(A35)}& {\displaystyle \frac{p{\mathrm{CO}}_{\mathrm{2},\mathrm{\infty }}}{p{\mathrm{CO}}_{\mathrm{2},\mathrm{0}}}}\approx \mathrm{1}+\mathrm{2}\cdot {\displaystyle \frac{\mathrm{\Delta }{N}_{\mathrm{o}}}{N_{\mathrm{o},\mathrm{0}}}}.\end{array}$$

From Fig. 3a we see that on multimillennial timescales the majority of the pulse is taken up by the ocean. Further, CaCO₃ compensation and weathering–burial imbalances in the CaCO₃ cycle (terrestrial neutralization) add (or remove) additional carbon to (from) the ocean. We thus assume that $\mathrm{\Delta }{N}_{\mathrm{o}}\approx \mathrm{1.5}\cdot 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₂

$$\begin{array}{}\text{(A36)}& {\displaystyle \frac{\mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2}}}{p{\mathrm{CO}}_{\mathrm{2},\mathrm{0}}}}\approx \mathrm{3}\cdot {\displaystyle \frac{P}{N_{\mathrm{o},\mathrm{0}}}}.\end{array}$$

We recall the definition of the airborne fraction:

$$\begin{array}{}\text{(A37)}& \mathrm{IRF}\left({\mathrm{CO}}_{\mathrm{2},\mathrm{a}}\right)=\mathrm{2.12}\text{GtC}\mathrm{\mu }{\text{atm}}^{-\mathrm{1}}\cdot {\displaystyle \frac{\mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2},\mathrm{a}}}{P}}.\end{array}$$

Finally, with Eqs. (A36) and (A37), we express the IRF at the time when CaCO₃ compensation is completed by

$$\begin{array}{}\text{(A38)}& {\mathrm{IRF}}_{{\mathrm{CaCO}}_{\mathrm{3}}}\left({\mathrm{CO}}_{\mathrm{2},\mathrm{a}}\right)=\mathrm{3}\cdot \mathrm{2.1}\text{GtC}\mathrm{\mu }{\text{atm}}^{-\mathrm{1}}\cdot {\displaystyle \frac{p{\mathrm{CO}}_{\mathrm{2},\mathrm{0}}}{N_{\mathrm{o},\mathrm{0}}}}.\end{array}$$

The remaining atmospheric fraction is thus independent of the pulse size for $P\ll N_{\mathrm{o},\mathrm{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₂ 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 δ¹³C perturbation

Carbonate compensation has almost a negligible influence on the δ¹³C perturbation. The carbonate compensation process acts, aside from small isotopic discrimination, equally on the different isotopes. Therefore, ¹²CO_2,a and ¹³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₃ adds an isotopic perturbation to the ocean–atmosphere–land system. The amount of carbon added corresponds to about half of the pulse size (P). CaCO₃ has a δ¹³C signature of about +2.9 ‰ compared to −24 ‰ of the pulse. We reference the CaCO₃ 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₃ dissolution is about 20 times smaller than that of the pulse $(\mathrm{0.5}\times P\times (\mathrm{2.9}\mathrm{‰}-\mathrm{0.8}\mathrm{‰}):P\times (-\mathrm{24}\mathrm{‰}+\mathrm{3}\mathrm{‰})=-\mathrm{1}:\mathrm{20})$. The dissolution of CaCO₃ enlarges the carbon reservoir by which the initial atmospheric δ¹³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₃ compensation on δ¹³C_a is small, and Eq. (A29) is still approximately valid for the timescale of CaCO₃ compensation.

Removal of the remaining CO₂ 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–ocean–land 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):

$$\begin{array}{}\text{(A39)}& \begin{array}{rl}& \mathrm{2}{\mathrm{CO}}_{\mathrm{2}\left(\mathrm{aq}\right)}+{\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\left(l\right)}+{\mathrm{CaSiO}}_{\mathrm{3},\left(\mathrm{s}\right)}\to {\mathrm{Ca}}_{\left(\mathrm{aq}\right)}^{\mathrm{2}+}\\ & +\mathrm{2}{\mathrm{HCO}}_{\mathrm{3}\left(\mathrm{aq}\right)}^{-}+{\mathrm{SiO}}_{\mathrm{2}\left(\mathrm{aq}\right)}.\end{array}\end{array}$$

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 δ¹³C perturbation by burial of organic matter and CaCO₃

In the case of δ¹³C, silicate rock weathering has no direct effect on the δ¹³C signature. The perturbation is removed on these long timescales by the replacement of carbon through burial of organic matter and CaCO₃ and weathering inputs:

$$\begin{array}{}\text{(A40)}& {\displaystyle \frac{d}{\mathrm{d}t}}(N_S\cdot \mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_S)=-F_{\mathrm{burial}}\cdot \mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{burial}},\end{array}$$

with N_S and ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\text{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 $\mathit{\tau }=N/F_{\text{burial}}(\sim \mathrm{90}$ kyr) (burial or weathering input is 0.46 GtC yr⁻¹; total amount: ∼40 000 GtC). Surprisingly, in the factorial simulation with sediments and constant weathering fluxes (SED₅₀₀), 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 δ¹³C perturbation is sensitive to Δδ¹³C gradients in the ocean. The isotopic signature of the burial flux is mainly determined by the upper-ocean signature. Changes in the gross exchange fluxes of carbon between the atmosphere and the ocean in response to perturbed CO₂ and changes in the isotopic fractionation of air–sea fluxes in response to perturbed temperatures increase the Δδ¹³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 Δδ¹³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 δ¹³C perturbation even on these very long timescales.

Appendix B: Sediment module and weathering feedback

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₃ and CaSiO₃ is parameterized as functions of temperature (T), runoff (R), and productivity (P) in the following form for CaCO₃:

$$\begin{array}{}\text{(B1)}& F_{{\mathrm{CaCO}}_{\mathrm{3}}}=F_{{\mathrm{CaCO}}_{\mathrm{3}},\mathrm{0}}\left(\mathrm{1}+k_{\mathrm{Ca}}(T-T_{\mathrm{0}})\right){\displaystyle \frac{R}{R_{\mathrm{0}}}}{\displaystyle \frac{P}{P_{\mathrm{0}}}}\end{array}$$

and for CaSiO₃

$$\begin{array}{}\text{(B2)}& F_{{\mathrm{CaSiO}}_{\mathrm{3}}}=F_{{\mathrm{CaSiO}}_{\mathrm{3}},\mathrm{0}}\cdot e^{\frac{\mathrm{1000}{E}_a}{R{T}_{\mathrm{0}}^{\mathrm{2}}}(T-T_{\mathrm{0}})}\cdot {\left({\displaystyle \frac{R}{R_{\mathrm{0}}}}\right)}^{\mathit{\beta }}{\displaystyle \frac{P}{P_{\mathrm{0}}}}.\end{array}$$

Subscript “0” indicates the initial value diagnosed prior to the pulse; k_Ca stems from a correlation between temperature and ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ ion concentration of groundwater with a value of 0.049; E_a is the activation energy for silicate weathering (63 kJ mol⁻¹); β is taken as 0.65; and R indicates runoff and is parameterized as a function of temperature:

$$\begin{array}{}\text{(B3)}& {\displaystyle \frac{R}{R_{\mathrm{0}}}}=\max \mathit{\{}\mathrm{0},\mathrm{1}+k_{\text{run}}(T-T_{\mathrm{0}}\left)\mathit{\right\}},\end{array}$$

with k_run=0.025. In the case of experiments 4box₅₀₀ and LGM₅₀₀, LPX-Bern was not coupled and instead productivity was parameterized, again following equations given in Colbourn et al. (2013):

$$\begin{array}{}\text{(B4)}& {\displaystyle \frac{P}{P_{\mathrm{0}}}}={\left({\displaystyle \frac{\mathrm{2}\frac{C}{C_{\mathrm{0}}}}{\mathrm{1}+\frac{C}{C_{\mathrm{0}}}}}\right)}^{\mathrm{0.4}},\end{array}$$

with C being atmospheric CO₂ 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).

Roth et al. (2014)Heinze et al. (1999)Heinze et al. (1999)Heinze et al. (1999)Heinze et al. (1999)Heinze et al. (1999)Roth et al. (2014)Jeltsch-Thömmes et al. (2019)Tschumi et al. (2011)Roth et al. (2014)Heinze et al. (1999)Jeltsch-Thömmes et al. (2019)Heinze et al. (1999)Anderson and Sarmiento (1994)Anderson and Sarmiento (1994)Anderson and Sarmiento (1994)Anderson and Sarmiento (1994)Roth et al. (2014)

Table B1 Parameters used for the sediment module, similar to what is shown in Appendix A of Tschumi et al. (2011) and with updates from Roth et al. (2014) and Jeltsch-Thömmes et al. (2019).