Negative carbon isotope excursions: an interpretive framework

Numerous negative carbon isotope excursions (nCIEs) in the geologic record occurring over 104–105 years are interpreted as episodes of massive carbon release. nCIEs help to illuminate the connection between past carbon cycling and climate variability. Theoretically, the size of a nCIE can be used to determine the mass of carbon released, provided that the carbon source is known or other environmental changes such as temperature or ocean pH can be constrained. A simple isotopic mass balance equation often serves as a first order estimate for the mass of carbon input, but this approach ignores the effects of negative carbon cycle-climate feedbacks. Here we show, using 432 earth system model simulations, that the mass of carbon release and associated environmental impacts for a nCIE of a given size and carbon source depend on the onset duration of that nCIE: the longer the nCIE onset duration, the greater the required carbon input in order to counterbalance the input of 13C-enriched carbon through carbonate compensation and weathering feedbacks. On timescales >103 years, these feedbacks remove carbon from the atmosphere so that the relative rise in atmospheric CO2 decreases with the nCIE onset duration. Consequently, the impacts on global temperature, surface ocean pH and saturation state are reduced if the nCIE has a long onset duration. The framework provided here demonstrates how constraints on the total nCIE duration and relative shape—together determining the onset duration—affect the interpretation of sedimentary nCIEs. Finally, we evaluate selected well-studied nCIEs, including the Eocene Thermal Maximum 2 (∼54 Ma), the Paleocene–Eocene Thermal Maximum (∼56 Ma), and the Aptian Oceanic Anoxic Event (∼120 Ma), in the context of our model-based framework and show how modeled environmental changes can be used to narrow down the most likely carbon emissions scenarios.


Introduction
Transient negative carbon isotope excursions (nCIEs) in the geologic record have garnered increasing attention because of the parallels with the modern climate experiment. As anthropogenic emissions of fossil fuel carbon accumulate in the atmosphere, the carbon isotopic composition of atmospheric carbon dioxide (δ 13 C CO2 ) has declined by ∼1.5‰ (the Suess effect) (Suess 1955, Keeling et al 1979, Keeling et al 1980, Mook 1986). Ultimately, the Suess effect will be recorded as a nCIE in the sedimentological record of the future (Norris et al 2013). Numerous nCIEs in the geological record are coeval with evidence for increasing global temperature and/or ocean acidification-a combination of impacts that are the hallmark of atmospheric greenhouse gas loading (Hönisch et al 2012). Investigating rates and masses of carbon release for past nCIEs in combination with environmental impacts informs us about climate sensitivity under various carbon emission scenarios and background climate states. Comparison of these events with the modern provides essential insight into the long-term response of the Earth system to fossil fuel emissions. An increasing number of nCIEs have emerged from the geological record as high-resolution measurements on marine sediment cores become more abundant, leading to hypotheses that these events may represent threshold or tipping point processes rather than independent, stand-alone events (Lunt et al 2011, DeConto et al 2012, Armstrong McKay and Lenton 2018. Hence, geological nCIEs may also provide an opportunity to investigate the sensitivity and role of carbon cycle feedbacks in amplifying or diminishing the consequences of massive carbon release. The carbon stored in the exchangeable reservoir, consisting of atmospheric carbon (CO 2 ), dissolved inorganic (DIC) and organic carbon (DOC) in the oceans, and carbon stored in land biota and soils, has a residence time of ∼10 5 years (figure 1). For nCIEs that are represented globally in both deep and shallow marine and terrestrial settings and occur on timescales shorter than 10 5 years, the presumption is these originate from the release of isotopically light carbon to the exchangeable carbon reservoir (Dickens et al 1995). Fluxes of carbon required to generate a given nCIE are often calculated by a simple isotopic mass balance approach (equation (1)). The mass of carbon required (M a ) to generate a nCIE of a particular magnitude (δ f -δ i ) in the exchangeable carbon reservoir (M i ) with an initial carbon isotopic composition (δ i ) depends on the source of the added carbon and hence its isotopic signature (δ a ). Volcanic activity, for instance, emits CO 2 with a δ 13 C that approximates the mantle value of −5 to −8‰ (Javoy et al 1986), while photosynthetic fractionation leaves marine and terrestrial organic matter more depleted in 13 C, with δ 13 C signatures ranging from −20 to −30‰ (Boutton 1991). Organic carbon reservoirs occurring in the form of coal, oil or gas may reach isotopic values lower than −40‰ (Boutton 1991) while microbial breakdown of organic matter drives the δ 13 C of the resulting methane even lower (more negative than −50‰) (Kvenvolden 1993). A particular-sized nCIE can be generated by the input of either a relatively large mass of carbon with higher δ 13 C or a small mass of carbon with lower δ 13 C. Hence, a well-known issue interpreting nCIEs is that the carbon isotopic record alone does not provide sufficient information to estimate the mass of carbon input responsible or the magnitude of change in atmospheric CO 2 associated with a nCIE. A greater mass of carbon added to the exchangeable reservoir should equate to a greater rise in atmospheric CO 2 on 10 3 -year timescales, with different implications for radiative forcing and temperature response.
A less often considered difficulty with the isotopic mass balance approach is that it assumes that the total mass of released carbon is instantaneously distributed through the exchangeable reservoir to produce a globally uniform nCIE. This assumption is problematic both for very rapid or very slow release of carbon. If carbon is released more rapidly than the mixing times between the various pools in the exchangeable reservoir (10 3 years), some of these carbon pools will record an amplified nCIE (e.g. Kirtland Turner and Ridgwell 2016). For The cGENIE exchangeable carbon reservoir in all our experiments includes 1770 Pg C atmospheric CO 2 with δ 13 C CO2 of −4.9‰, 30 300 Pg C DIC with δ 13 C DIC of 1.7‰, 7 Pg C of marine dissolved organic carbon with a δ 13 C of −27.4‰. Note that the version of cGENIE used here does not include a representation of terrestrial biomass or soil carbon and there is no explicit standing stock of marine biomass. carbon release over longer timescales that approach the residence time of carbon in the exchangeable reservoir, it becomes crucial to account for feedbacks such as carbonate compensation and rock weathering. More recently, modeling approaches have been utilized that explicitly consider the timing of a nCIE derived from sedimentary age models (Cui et al 2011, Cui et al 2013, Gutjahr et al 2017, Dunkley Jones et al 2018 in calculating the mass of carbon required to generate a particular nCIE. However, even constraining the size of a nCIE to use in isotope mass balance calculations or more sophisticated modeling for a single event is not straightforward. Marine carbonates are often used as targets for reconstructing nCIEs because they are presumed to reflect the isotopic composition of the DIC pool, by far the largest pool in the exchangeable surface reservoir (figure 1). However, marine carbonate nCIEs may vary between locations for the same event due to many factors, including spatial variations in bioturbation and dissolution intensity (Sluijs and Dickens 2012) and whether the carbonate originated from a benthic (seafloor) or planktonic (surface) source. Water depth (for benthic carbonates) may also lead to variations between records. It is possible to identify nCIEs from organic matter as well, although these records must be interpreted with the further caveat that changes in atmospheric CO 2 may lead to changes in the fractionation between DIC and organic matter or changes in organic matter composition that overprint the record of the nCIE as seen by the global exchangeable reservoir (Corsetti et al 2005).
When considering multiple nCIEs in the geologic record, there is wide variation not only in magnitude, but also in shape and duration, all of which complicate the calculation of required carbon input. This nuance is not captured in a simple isotopic mass balance approach. Importantly, the goal in using isotopic mass balance to infer the magnitude of carbon released is often to constrain past climate sensitivity (or the global temperature increase due to a doubling of atmospheric CO 2 ) by relating this carbon release to a change in atmospheric CO 2 and combining this with estimates of temperature change. Uncertainty in carbon source and nCIE size will result in the misinterpretation of carbon input required to explain δ 13 C excursions (equation (1)) but assumed duration and relative timing of an event will also influence (1) the carbon input required, as various reservoir exchanges and feedbacks act on different timescales and (2) the proportional increase in atmospheric CO 2 by determining the fraction of the total emitted carbon that remains in the atmosphere following its release. Hence, uncertainty in age models translates into uncertainty in reconstructed climate sensitivity because for a given mass of carbon release to the exchangeable reservoir, one cannot infer the increase in atmospheric CO 2 without constraining the timing of that carbon input. Instantaneous carbon release will lead to the maximum change in atmospheric CO 2 while slower release of the same mass of carbon will dampen the atmospheric CO 2 rise.
The many nCIEs across the geological record demonstrate that the carbon cycle has been disrupted frequently throughout Earth's history (Saltzman and Thomas 2012). Perhaps the most well studied nCIE, the Paleocene-Eocene Thermal Maximum (PETM, ∼56 Ma) (figure 2(b)), is characterized by a rapid onset of several thousand years, a negative δ 13 C excursion of ∼3-4‰ and complete isotopic recovery on timescales of 10 5 years (Röhl et al 2007, Zachos et al 2008, McInerney and Wing 2011, Westerhold et al 2018a. Another well-known nCIE is the Eocene Thermal Maximum 2 (ETM-2) (figure 2(c)), which is characterized by a nCIE that is half the magnitude of the PETM and more symmetrical in shape, i.e. the onset duration is equal to the duration of the recovery phase (Lourens et al 2005, Stap et al 2009. The PETM, ETM-2, and other smaller nCIEs of the early Cenozoic are known collectively as hyperthermals (Thomas and . The red lines indicate simulated nCIEs from this study that most closely resemble the bulk δ 13 C records. Zachos 1999). Many of these purported hyperthermals are more similar to the magnitude and shape of the ETM-2, i.e. nCIEs of 1-2‰ occurring over 40-100 kyr (Sexton et al 2011, Kirtland Turner et al 2014, Westerhold et al 2018b. Though a precise definition of 'hyperthermal' is lacking, these events are typified by simultaneous deep-sea warming and carbonate dissolution in addition to a negative δ 13 C excursion. While not always described in the 'hyperthermal' terminology, a number of Mesozoic events share similar characteristics. For instance, the Toarcian Oceanic Anoxic Event (T-OAE, ∼183 Ma) is characterized by a nCIE of ∼6‰ between two positive δ 13 C excursions (Jenkyns 1988, Hermoso et al 2012, Müller et al 2017. Similarly, nCIEs have been reconstructed in association with the end-Permian mass extinction (∼252 Ma) (Meyer et al 2011), the end-Triassic mass extinction (∼201 Ma) (Pálfy et al 2001), and preceding multiple Cretaceous OAEs (particularly Aptian OAE-1a (∼120 Ma) (Menegatti et al 1998) (figure 2(a)), and OAE-1b (∼111 Ma) (Wilson and Norris 2001). One complication in comparing these older events to the Cenozoic hyperthermals is that deep-sea records are scarce, and nCIEs are often reconstructed from shallow marine settings using carbonate that may have an unclear diagenetic history or using organic matter that may reflect a combination of sources (e.g. Corsetti et al 2005). However, carbon isotope records of these events have been previously evaluated as reflecting changes in the exchangeable carbon reservoir as a result of massive carbon input similar to the Cenozoic hyperthermals.
Given the range of nCIEs observed in the geologic record in terms of size, duration, and shape, as well as varying confidence regarding age control, we herein present an ensemble of modeled nCIEs using the Earth system model 'cGENIE' to reproduce variable nCIEs assuming various carbon sources. We generate 432 experiments (figure 3) and calculate the associated fluxes of carbon and impacts on marine and atmospheric chemistry and temperature for each. Advancing beyond a simple isotope mass balance approach, our modeling utilizes a 3D dynamical ocean model and incorporates carbon cycle feedbacks such as changes in ocean solubility and carbon speciation, carbonate compensation, and temperature-dependent rock weathering. Our framework facilitates assessment of individual nCIEs throughout the geologic record by providing an interpretive template for relating nCIEs of a given size, duration, and shape in the marine DIC reservoir to particular carbon cycle drivers and climatic and environmental consequences.

Methods cGENIE model and set-up
In this study, we use the intermediate complexity carbon-Grid ENabled Integrated Earth system model (cGENIE). cGENIE comprises a 3D frictional geostrophic ocean model and a simple energy-moisture balance atmosphere (Edwards and Marsh 2005). We include a CO 2 -climate feedback: radiative forcing due to CO 2 doubling in cGENIE is ∼4 W m −2 . We also include temperature-dependent continental weathering (Colbourn et al 2013) with a 1:1 carbonate to silicate weathering ratio and allow transport of alkaline species from the terrestrial into the marine environment. To balance terrestrial input, we couple a sediment model that calculates the sedimentary preservation of inorganic carbon in the form of calcium carbonate (CaCO 3 ) in deep sea sediments . Biogeochemical cycling of elements and isotopes is represented throughout these four model components . A detailed description and evaluation of the carbon isotope (δ 13 C) cycle in cGENIE is provided in  and Kirtland Turner and Ridgwell (2016). Marine (export) productivity is calculated by the local availability of phosphate with modifiers for light limitation and sea ice extent . Phosphate is modeled as a closed system, with no input or burial fluxes, so its oceanic distribution is controlled by a combination of modeled patterns in export production and physical ocean circulation (which returns nutrients from the deep ocean to the surface). Calcium carbonate (CaCO 3 ) export is calculated using a fixed ratio of calcium carbonate to particulate organic carbon (POC) of 0.2. Both POC and CaCO 3 are remineralized in the water column using fixed depth-dependent profiles. Any POC reaching the seafloor is remineralized there (hence we do not account for sedimentary organic carbon burial). Importantly, our model also excludes any representation of land biota or soils, so that the exchangeable carbon reservoir is approximated as atmospheric carbon and marine dissolved inorganic and organic carbon only, which make up about 95% of the preindustrial exchangeable reservoir (figure 1).
Because our goal is to provide a general interpretative framework for nCIEs applicable throughout the geological record, we use a simplified low-resolution (18×18) continental configuration with a single symmetric pole-to-pole continent and 16 ocean depth levels. Still, the model ocean circulation is comparable to that in a higher resolution configuration run for the late Paleocene (figure S1, tables S1-S2 is available online at stacks.iop.org/ERL/14/085014/mmedia). Both configurations show a dominant Southern Ocean overturning circulation, though the configuration used here shows overall weaker overturning strength and more symmetrical circulation between hemispheres. Due to differences in the overturning circulation and hence differences in nutrient distribution, the symmetric configuration yields slight increases in fluxes of particulate organic and inorganic carbon and a slight increase in the size of the marine dissolved organic carbon reservoir. However, due to the highly idealized continental configuration and its influence on the modeled circulation pattern, we cannot make explicit comparisons to the behavior of ocean circulation in previous Paleocene cGENIE experiments (e.g. Kirtland Turner and Ridgwell 2016).
The model is equilibrated to greenhouse climate conditions, with atmospheric CO 2 concentrations of three times preindustrial values (834 ppm), to approximate Mesozoic and Cenozoic greenhouse conditions (Fletcher et al 2008). We use the same initial carbon cycle conditions as the published cGENIE Paleocene configuration including: atmospheric δ 13 C of −4.9‰, adjusted major ion concentrations (Ca 2+ 18.2 mmol kg −1 , Mg 2+ 29.9 mmol kg −1 ), and lower alkalinity (1975 μmol kg −1 ) along with modern (preindustrial) inventories for phosphate (2.16 μmol kg −1 ) and sulfate (15.0 mmol kg −1 ) (Ridgwell and Schmidt 2010). The model is run to steady state with carbon input from volcanic outgassing and alkaline run-off from the continent balanced by carbon removal through burial of CaCO 3 . The resulting modeled global average sea surface temperature is ∼26°C, sea ice is absent and mean sedimentary wt% CaCO 3 is 42%.
We employ an inverse modeling technique to produce nCIEs wherein the isotopic signature of the surface ocean dissolved inorganic carbon pool (δ 13 C DIC) is forced to follow a prescribed δ 13 C evolution, while an internal algorithm determines how much atmospheric CO 2 input with a specified isotopic signature is necessary to match this profile (Cui et al 2011, Cui et al 2013, Kirtland Turner and Ridgwell 2013, Dunkley Jones et al 2018. We choose to force surface DIC δ 13 C rather than atmospheric CO 2 δ 13 C because the former is a closer approximation to data from marine carbonates. When modeled surface DIC δ 13 C values are higher than prescribed, a pulse of isotopically light CO 2 is released into the atmosphere. CO 2 with the same isotopic signature is removed from the atmosphere when modeled surface DIC δ 13 C drops below the prescribed value. Hence, CO 2 is artificially removed when carbonate compensation and weathering feedbacks are unable to sufficiently restore surface DIC δ 13 C values, but this carbon removal does not represent any particular physical process. We impose carbon input (and removal) with isotopic signatures of −6‰, −12‰, −22‰ and −60‰, conceptually corresponding to carbon input from volcanism, a mixture of volcanism and organic carbon, organic carbon, and biogenic methane, respectively. The modeled nCIEs vary in size (0.5‰, 1.0‰, 3.0‰ and 6.0‰) and duration (50, 100 and 300 kyr) to represent the wide range of nCIEs recorded in the geological record. Furthermore, we define nine nCIE shapes (figure 3). We identify three dimensions of symmetry in duration: events can be (1) symmetrical in duration with equal duration onset and recovery phases (δ 13 C minimum occurs halfway through the event duration) or asymmetrical with either (2) a relatively rapid onset and slow recovery (δ 13 C minimum occurs at ¼ of event duration) or (3) relatively slow onset and rapid recovery (δ 13 C minimum occurs at ¾ of event duration). We also identify three dimensions of symmetry in size: (1) δ 13 C values at the onset and recovery are identical, (2) final δ 13 C values are 50% higher than initial (overshoot) and (3) final δ 13 C values are 50% lower than initial (undershoot). Hence, the most rapid modeled nCIE onset is 12.5 kyr in experiments with a total duration of 50 kyr. The slowest nCIE onset of 225 kyr occurs in experiments with a total duration of 300 kyr. A 1‰ nCIE with an overshoot will end with DIC δ 13 C values 0.5‰ higher than initial while a 1‰ nCIE with an undershoot will end with DIC δ 13 C values 0.5‰ lower than initial.

Results
Our large ensemble of nCIE experiments allows us to compare both the required carbon input and removal in terms of total mass and fluxes for particular nCIEs and relate the range of these values to uncertainty in carbon source, event duration, and event shape. We calculate 1000-year-binned average values of input and removal fluxes and sum net masses of carbon input/removal over the onset/recovery phases. We also compare the resulting carbon cycle and climate changes in the context of these same nCIE characteristics. In particular, we evaluate modeled changes in atmospheric CO 2 , sea surface temperature (SST), surface ocean pH, surface ocean saturation state, weathering rates and wt% CaCO 3 at the nCIE peak.
Diagnosed carbon forcing All modeled nCIEs require an input of isotopically light carbon during the onset phase, while carbon needs to be removed during nCIE recovery to restore surface DIC δ 13 C to higher values. The total mass of carbon necessary to produce a nCIE of specific size depends to a first order on the isotopic signature of the source carbon. Total masses of carbon input diagnosed in our experiments vary by three orders of magnitude (figures 4, S2). On the high end, to generate the largest size (6‰) nCIE from volcanic carbon with a δ 13 C of −6‰ requires a carbon input mass that is an order of magnitude larger than the entire mass of ∼32 000 Pg C in the modeled exchangeable reservoir at the start of the experiment. This contrasts with just a few thousand Pg C necessary to generate the same nCIE with carbon derived from methane with δ 13 C signature of −60‰ (figures 4(a), S2). Significant variability in the required carbon input is also generated by variations in the duration and shape of the nCIE alone. Our results generally show that increasing assumed nCIE duration leads to an increase in the total carbon input on the 10 4 -year input timescales modeled here. Increases greater than 50% occur for some combinations of nCIE size and carbon source (figure 4(a)-left panel and figure S2-moving left to right within each of subplots e-h). Assuming the total duration of the nCIE is also known, varying the duration symmetry still causes variability in the calculated total carbon input. Increasing the relative onset duration from rapid to slow for a nCIE of a given size and total duration generally increases the total carbon input required (figure 4-compare 4(b) to 4(c), figure S2-moving from the top to bottom row within each column, i.e. comparing subplots a, e, and i for a given nCIE size, figure S3). The effect is more pronounced for longer total durations. Whether a nCIE ends with an undershoot or an overshoot in surface DIC δ 13 C does not impact the total mass of carbon input required during the onset phase.
Total carbon removal necessary to restore surface DIC δ 13 C for a nCIE of a given size and shape is similarly controlled to a first order by the imposed isotope composition of the removed carbon. In all our experiments, additional carbon removal occurs via the same isotopic composition as carbon input, allowing us to make a straightforward comparison between the relative masses of carbon input and removal, which we illustrate as the net carbon input to the exchangeable reservoir ( figure S4). nCIEs with a δ 13 C overshoot (figure S4, left columns) are generally created by net removal of carbon over the total nCIE duration (indicated by blue colors) except for those with a long total duration of 300 kyr. nCIEs ending with a δ 13 C undershoot (figure S4, right columns), or nCIEs that are symmetrical in size ('no-shoot', figure S4, middle columns) always require a larger mass of carbon input compared to carbon removal (indicated by red colors). In other words, the average carbon isotopic composition of the surface DIC pool is exactly restored, but to achieve this, less carbon is removed than is added. Similarly, nCIEs that are symmetrical in size or with an undershoot leave the exchangeable reservoir significantly larger than at the experiment onset, particularly in the case of the largest nCIEs with the longest duration. In contrast, nCIEs with an overshoot always leave the exchangeable reservoir smaller than at the experiment onset, regardless of nCIE size or duration (figure S5).
Rate of carbon input to generate a nCIE of a given size also depends to a first order on the source of carbon input (figures 4, S5). The maximum sustained rates of carbon input, calculated from 1000-year bins of the yearly average rate, exceed 20 Pg C yr −1 and are required to generate the 6‰ nCIE using volcanic carbon with δ 13 C of −6.0‰. Total event duration is significant for determining rate of carbon input for a nCIE of a given size and carbon source. Decreasing total duration from 300 to 50 kyr more than doubles the maximum rate of carbon input for all nCIE sizes (figure S6-moving from right to left within each subplot). In other words, compared with the total carbon input, the rate of carbon input is more sensitive to calculated event duration. Rate is also sensitive to assumptions about duration symmetry. For a nCIE of a given duration, varying the shape from a slow onset to a rapid onset can result in more than a doubling of the maximum sustained rate of carbon emissions. Hence, the relationship between nCIE duration or shape and maximum carbon flux is the opposite as between nCIE duration or shape and gross carbon input ( figure 4(a)). Maximum carbon flux is highest for shorter, rapid onset events whereas gross carbon input is highest for longer, slower onset events. Overall, our ensemble results in an enormous range of maximum sustained rates in carbon input during nCIE onset phases. The smallest nCIEs require rates of ∼0.01 Pg C yr −1 assuming depleted carbon sources of either biogenic methane or organic carbon. However, even the smaller nCIEs of 0.5 and 1.0‰ require relatively large sustained rates of carbon input (between 0.03 and 1.91 Pg C yr −1 ) if the assumed carbon source is volcanism with a δ 13 C of −6‰. None of our simulated nCIEs forced with biogenic methane (δ 13 C −60‰) or organic carbon (δ 13 C −22‰) require rates of carbon input reaching the current anthropogenic emission rate of 10 Pg C yr −1 . Even a 6‰ nCIE occurring over 50 kyr with a rapid onset over just 12.5 kyr requires a maximum sustained input of only about 1.4 Pg C yr −1 assuming an organic carbon (δ 13 C −22‰) source. In fact, only the experiments with volcanic (δ 13 C −6‰) carbon as the source ever require input rates in excess of the current anthropogenic value. All nCIEs forced with the mixture of volcanism and organic carbon (δ 13 C −12‰) have maximum diagnosed emissions less than 6 Pg C yr −1 .

Carbon cycle and climate impacts
We focus on carbon cycle and climate changes that could feasibly be estimated using available proxies, including the change in atmospheric CO 2 (figures 5, S7), surface ocean pH (figure S8), surface ocean saturation state (figure S8) and SST (figure S9). Given the dynamic ocean model in cGENIE, ocean circulation in our experiments is not fixed, yet we find minor changes of less than ∼3% in ocean overturning at the nCIE peak for most experiments (less than 1 Sv for all experiments except the 6‰ nCIEs forced with −6‰ volcanic carbon, which still show less than 3 Sv change at the nCIE peak). These subtle changes in ocean circulation also result in small changes in export production, which, in our model, is a function only of phosphate availability with modifiers for light limitation and sea ice extent (although no sea ice occurs in our experiments). These export production changes are small relative to changes in the size of the atmospheric and inorganic ocean carbon reservoirs (less than 1 Pg C yr −1 change at the nCIE peak in all but a few 6‰ nCIE experiments that show a maximum change of 3 Pg C yr −1 at the nCIE peak). The simulated changes to the carbon cycle are thus driven primarily by the total mass and rate of carbon input and associated feedbacks in the carbonate chemistry-the focus of our discussion henceforth-rather than changes in primary productivity or ocean circulation.
Our experiments result in an extremely large range of rates and total masses of carbon input (figure 4) that diverge from estimates based on the simple mass balance approach (figure 5). The most pronounced difference between the mass balance approach and our nCIE simulations occurs for the nCIE of 6‰ with a 225 kyr-long onset phase forced with δ 13 C carbon of −6‰. Simple isotope mass balance predicts that ∼147 200 Pg C input is required, while cGENIE requires nearly three times as much carbon input (∼413 800 Pg C). We use the results for environmental variables, particularly temperature, to test whether these results are realistic. For instance, the maximum atmospheric CO 2 reached in any of our experiments is over 120 000 ppm CO 2 , which is equivalent to slightly more than seven doublings from the initial concentration of 834 ppm. In cGENIE, the imposed radiative forcing due to CO 2 doubling is ∼4 W m −2 resulting in a typical surface temperature rise of 3.2°C per CO 2 doubling (figure S12). The largest increases in atmospheric CO 2 occur in our long duration (300 kyr) experiments with a 6‰ nCIE driven by a −6‰ (volcanic) source, leading to an increase in average surface air temperature of almost 28°C or an absolute maximum average surface air temperature of ∼50°C. SST rises by ∼22°C. In contrast, when the most isotopically depleted carbon source is assumed, this yields an increase in atmospheric CO 2 equivalent to approximately one doubling and SST warming of just 2.6°C. The smallest nCIEs of 0.5‰ lead to a minor increase in atmospheric CO 2 of less than 150 ppm if the assumed carbon source is either organic carbon or biogenic methane. The associated increase in SST is just a few tenths of a degree Celsius. For a nCIE of a given size and carbon source, the maximum increase in CO 2 and maximum increase in SST both correspond to the maximum carbon addition rate rather than the total mass of carbon added, such that shorter durations and rapid onsets lead to greater increases in CO 2 and temperature. The sole exception occurs for 6‰ nCIEs forced by volcanic carbon (−6‰), where larger increases in CO 2 and SST occur for relatively longer duration events with slower onsets. For these nCIEs only, the maximum change in CO 2 and SST corresponds to the largest gross carbon input rather than the highest carbon input rate. The atmosphere in these scenarios becomes a proportionally larger carbon reservoir than the ocean, i.e. the ocean contains 17 121 Pg C at the nCIE peak, while the atmosphere contains nearly double that mass (27 429 Pg C) as a consequence of reduced CO 2 solubility with extreme warming. This means that ocean uptake is unable to remove excess CO 2 even when that CO 2 is emitted relatively slowly, so it continues to accumulate in the atmosphere. It is important to note that such large masses and rates of carbon input may not be representative of the events that occurred during any time in the Mesozoic or Cenozoic. With greater atmospheric carbon input, greater quantities of CO 2 dissolve into the surface ocean, lowering surface ocean pH and the saturation state with respect to calcite and aragonite (Ω cal and Ω arg , respectively). The largest change in mean surface ocean pH (a decrease of 1.7) occurs in 6‰ nCIEs forced with volcanic carbon. Aside from the 6‰ nCIEs forced by enormous quantities of volcanic carbon described above, the relative decline in surface ocean pH for each experiment is equivalent to the relative increase in CO 2 -in other words, a shorter onset duration and hence more rapid rate of carbon input leads to a larger decline in surface ocean pH for a nCIE of a given size and carbon source. The largest change in saturation state with respect to calcite (a decrease of 5) occurs in 6‰ nCIEs forced with volcanic carbon that have the shortest onset durations of 12.5-25 kyr. In the initial model state prior to carbon emissions, the surface ocean is everywhere oversaturated with respect to calcite (Ω cal >1) and 99.2% of the grid cells are oversaturated with respect to aragonite (Ω arg >1). While aragonite undersaturation (Ω arg <1) can occur in parts of the surface ocean with very modest carbon forcing, calcite undersaturation requires input greater than 10 000 Pg C. Both calcite and aragonite undersaturation occurs predominately at the high latitudes where carbon uptake is maximized due to the higher solubility of CO 2 . Carbon forcing greater than 40 000 Pg C over 12.5 kyr or shorter is required to drive Ω cal <1 at low latitudes.
Onset duration (or the time to the nCIE peak) controls whether or not the response of the mean surface ocean pH and saturation state behave similarly in response to a given carbon input (figure S8). As nCIE duration increases, the magnitude of the surface pH decrease is slightly reduced (with the exception of 6‰ nCIEs forced by volcanic carbon (−6‰)), but the magnitude of the surface saturation state decrease is more notably reduced. The same relative effect is seen as a function of temporal symmetry-a rapid onset and slow onset of a nCIE with the same total duration show similar declines in surface ocean pH but the rapid onset nCIE shows a larger decline in saturation state. Hence, the mean surface ocean pH is more sensitive to total C input and surface saturation state is more sensitive to C flux. These results are consistent with the findings of Hönisch et al (2012) who demonstrated that ocean pH and saturation state change are progressively decoupled when carbon input occurs over long timescales. Notable decoupling begins for carbon input over 10 4 years or more.
We can also evaluate how each of the above variables evolves with time across each experiment. As an example, we show the response of atmospheric CO 2 , surface ocean pH and saturation state with respect to calcite, deep sea carbonate content (wt% CaCO 3 ), weathering flux (in terms of Ca 2+ ) and SST for a modeled nCIE similar to the Eocene Thermal Maximum 2 (ETM-2, 53 Ma) ( figure 6(a)). We take a nCIE of 1‰ with a total duration of 100 kyr and a symmetric shape (50 kyr onset duration) as the best approximation and evaluate the modeled response of the above variables across the experiment duration. Atmospheric CO 2 and SST peak at the same time as gross C input, coincident with minimum δ 13 C in surface ocean DIC. Simultaneously, surface ocean pH and saturation state reach their minima while weathering fluxes reach their maxima. Sedimentary wt% CaCO 3 , however, shows a minimum that precedes the peak in SST and weathering flux and has already begun to rebound before the δ 13 C minimum. The later rebound in CaCO 3 shows maximum accumulation approximately coeval with maximum surface ocean saturation state. In a similar manner, a comparison is made between the PETM, OAE-1a, and the nCIE simulations performed to establish the framework presented here (figure S11).

Diagnosed carbon input
Our experiments highlight the significance of considering timing in diagnosing the carbon input required to generate a nCIE. Of particular importance is the time to peak nCIE (onset duration), which is controlled by both nCIE duration and nCIE shape in our experimental framework. Overall, we find that greater carbon input is necessary to generate a nCIE with a longer onset duration, but the rise in atmospheric CO 2 decreases with onset duration. The experiments start from steady state with respect to the exchangeable carbon reservoir (carbon inputs from weathering fluxes and volcanism are balanced by sedimentary carbonate burial), on top of which carbon inputs diagnosed by the inversion scheme are added from an external, specified source. Carbon addition diagnosed by the inversion scheme during the nCIE onset directly increases the size of the exchangeable carbon reservoir. The exchangeable reservoir also grows or shrinks because of the indirect effects of global temperature and ocean chemistry on weathering rates and the preservation of sedimentary CaCO 3 . These feedbacks operate over characteristic timescales, such that nCIE duration exerts an important control on the total carbon mass in the exchangeable reservoir as well as the relative distribution of carbon among atmospheric CO 2 and the oceanic DIC and DOC pools.
At the start of each experiment, the exchangeable reservoir totals just over 31 800 Pg C (atmospheric CO 2 of ∼1800 Pg C, marine DIC of ∼30 000 Pg C, and marine DOC of ∼6.7 Pg C). For a 6‰, 300 kyr nCIE with a relatively rapid onset (75 kyr) driven by −22‰ carbon, the masses of these reservoirs total just over 48 000 Pg C (atmosphere ∼7000 Pg C, DIC ∼41 000 Pg C, and DOC ∼6.9 Pg C) at the nCIE peak, an increase of about 17 000 Pg C over initial conditions. However, the inversion scheme diagnoses a gross C input over this interval of 13 760 Pg C. Thus >3000 Pg C (more than 20% of the total exchangeable carbon reservoir increase) is attributable to excess weathering not balanced by increased CaCO 3 burial by the time of the nCIE peak. Further, there has been a relative redistribution of carbon such that the atmosphere is a proportionally larger carbon pool (although DIC still dominates the exchangeable reservoir mass) because as atmospheric CO 2 rises and surface ocean pH declines, the ocean's buffering capacity is reduced.
A mass balance framework utilizes only information about the initial mass of the exchangeable carbon reservoir, the size of the nCIE, and the inferred carbon source. Our more nuanced framework, by incorporating inorganic carbon cycle feedbacks, also utilizes constraints provided by age models on the total and relative duration of a nCIE and provides a refined first order estimate regarding the range of possible carbon fluxes and cumulative input masses. Figure 5 illustrates the difference between the simple isotope mass balance calculations for carbon input and our cGENIE model results for gross carbon addition. While the simple isotope mass balance approach (equation (1)) provides a single value for a particular size nCIE given an assumed carbon source (with specified δ 13 C), our experiments generate a range of estimates for carbon input, depending on the duration and shape of the nCIE. The discrepancy between our cGENIE diagnosed carbon inputs and a simple mass balance approach is most pronounced for larger nCIEs generated using carbon sources with higher δ 13 C. For instance, assuming a volcanic carbon source (figure 5(a)), the isotope mass balance approach provides an estimate of input mass that falls nearly in the middle of the range calculated by cGENIE for a 0.5‰ and 1‰ nCIE, but isotope mass balance gives a value closer to the lower range of diagnosed carbon inputs for nCIEs of 3‰ and 6‰. This increasing discrepancy between isotope mass balance and our results for large input masses is due to the temperature-dependent weathering parameterization (i.e. weathering increases by more and thus gross carbon must increase by relatively more to compensate as total carbon input increases).
Both our approach and the simple isotope mass balance (equation (1)) require information about initial carbon cycle conditions. In the mass balance approach, this is simply the mass of the total exchangeable reservoir. Our approach requires not only this information, but also makes presumptions about the initial buffering capacity of the oceans, which will influence how carbon is exchanged between reservoirs and how much seafloor carbonate dissolves. These carbon cycle assumptions must be considered when comparing our results against particular geological events, especially those that occurred prior to the advent of widespread pelagic calcification during the Jurassic (∼170 Ma) (Arvidson et al 2014).
We do not explicitly test the impact of varying initial carbon cycle conditions on our results; however, we can infer the general impact of particular assumptions. While it is likely that the marine DIC pool (by far the largest of the exchangeable carbon pools) was similar in mass through much of the Cenozoic, it was likely larger during much of the Mesozoic (e.g. Ridgwell 2005). If the mass of the total exchangeable reservoir (atmospheric CO 2 and ocean DIC and DOC) was significantly larger, the required carbon input of a given isotopic signature for a given nCIE would have been larger. Varying initial carbon cycle conditions, including alkalinity or major ion composition, would have a more nuanced influence, altering the relative partitioning of carbon between reservoirs and controlling the strength of feedbacks that are sensitive to this distribution. For example, initial carbon cycle conditions characterized by higher relative partitioning into atmospheric CO 2 (lower surface ocean pH) requires relatively larger masses of CO 2 emitted to achieve the same nCIE in the surface ocean DIC pool because a greater proportion of emissions will remain in atmosphere, i.e. the global mean Revelle factor is higher in a high pCO 2 /low pH world. With a greater increase in atmospheric CO 2 , the resultant enhanced weathering that releases relatively isotopically heavy alkaline carbon species will dampen the signal of the 13 C-depleted emissions and hence require even greater 13 C-depleted CO 2 emissions to generate the same nCIE in the surface ocean. Thus, when interpreting an individual geologic event in the context of our hypothetical framework, it is necessary to consider whether initial mean surface ocean pH is likely to have been lower or higher than what is modeled here (∼7.7).

Diagnosed carbon removal
Restoration of surface DIC δ 13 C to pre-excursion values during the recovery phase requires the addition of relatively heavy carbon (via carbonate dissolution and weathering) and/or removal of isotopically light carbon. The model calculates fluxes of carbon (and its isotopes) both within the exchangeable reservoir and including sediment and weathering feedbacks, and then compensates for any remaining deviations of the surface ocean DIC δ 13 C pool from the target curve through removal of CO 2 with a specified δ 13 C from the atmosphere. In other words, carbon simply disappears from the model atmosphere if the simulated carbon cycle feedbacks have failed to restore surface ocean DIC δ 13 C sufficiently to match the prescribed curve. This process is conceptually simple, but it has consequences that are most likely not realistic. First, this diagnosed additional atmospheric carbon removal is not simulated as any particular physical mechanism. Second, the fact that the carbon removal flux has a specified isotopic composition (here we use the same value as the carbon input), means that we are not adjusting the isotopic composition of the removed carbon based on the composition of exchangeable reservoir or based on likely fractionations of real removal processes (e.g. organic matter formation). Third, our method does not account for removal of carbon at a different isotopic composition than the input, though in reality these values may be very different (e.g. input of volcanic CO 2 followed by removal via organic carbon burial). An alternative approach for calculating carbon removal, even without explicitly modeling organic carbon burial, would be to remove carbon with an isotopic composition matching that of the modeled particulate organic carbon pool (e.g. Gutjahr et al 2017), rather than using an assigned value.
Only nCIEs with an undershoot where the final δ 13 C is lower than the initial δ 13 C, and/or a long duration and slow recovery result in diagnosed carbon removal fluxes substantially less than input fluxes. In other words, they show a large positive net carbon input ( figure S4). This indicates that modeled feedbacks are nearly sufficient to match the surface DIC δ 13 C signature in these cases and little additional diagnosed carbon removal is necessary. In modeled nCIEs where the diagnosed mass of carbon removal rivals or even exceeds the diagnosed carbon input, this implies that the carbonate compensation and weathering feedbacks alone were insufficient to drive the recovery of surface DIC δ 13 C and additional carbon removal (modeled as CO 2 artificially removed from the atmosphere) was necessary. A deficiency in our experiments is thus the lack of organic carbon burial, which might help explain the recovery of surface ocean DIC δ 13 C without requiring an additional removal flux via our inversion methodology for the larger and/or shorter duration nCIEs. Given that a number of nCIEs in the geologic record are followed by black shale deposition and associated positive δ 13 C excursions, e.g. the Toarcian and Aptian OAEs, there is clear evidence that organic carbon feedbacks not only played a role in recovery of the exchangeable reservoir from carbon input but often led to δ 13 C 'overshoot' behavior (Jenkyns 2010).

Carbon cycle and climate impacts
The change in atmospheric CO 2 and associated warming generated in each of our experiments allows us to reasonably rule out certain carbon input scenarios for certain nCIEs. In particular, the massive carbon inputs required to generate 6‰ nCIEs using volcanic carbon of −6‰ raise atmospheric CO 2 to levels that lead to an average ocean temperature greater than 40°C, inconsistent with data for even the Mesozoic greenhouses. The highest reconstructed low latitude surface temperatures during the Cretaceous are reconstructed from the TEX 86 archaeal lipid SST proxy, and are a maximum of ∼36°C (Schouten et al 2003). Further, evidence for SST in excess of 40°C during the PETM from Tanzania in combination with the exclusion of planktic foraminifera has been provided as evidence of a temperature threshold for these organisms (Aze et al 2014). It seems likely that widespread temperatures this high would have resulted in significant extinction and widespread dead zones. Adding the effects of severe thermal stress on marine organisms (not modeled here) could increase CO 2 even further if the biological carbon pump were to collapse. This suggests that large nCIEs occurring over the timescales modeled here and similar to several Mesozoic nCIEs are unlikely to have been driven by isotopically heavy carbon input (e.g. mantle-derived volcanism) alone. nCIEs of 3‰ forced with −6.0‰ carbon and 6‰ nCIEs forced with −12‰ carbon have maximum CO 2 concentrations 7500-17 000 ppm and maximum global mean SSTs of 33°C-36°C, more consistent with proxy temperature records.
Our experiments highlight the importance of nCIE duration and shape for determining carbon input rate, which has significant impacts on atmospheric CO 2 , SST and saturation state (figure S10). While nCIEs with short onset durations require less total input of carbon than nCIEs with longer onset durations, maximum atmospheric CO 2 concentrations are highest in the shortest nCIEs. Again, this can be attributed to marine carbonate compensation and terrestrial rock weathering, which both act to stabilize ocean saturation state and modulate atmospheric CO 2 concentrations over 10 4 -10 5 -year timescales (Walker et al 1981, Archer and Maier-Reimer 1994, Archer et al 1997, Archer 2005, Archer et al 2009a, Goodwin and Ridgwell 2010. In particular, the temperature-dependent weathering feedback in cGENIE operates with an e-folding timescale of ∼240 kyr (Colbourn et al 2015). As carbon input extends over intervals of time commensurate with the e-folding weathering timescale, the required total carbon input increases while the effect of those emissions on atmospheric CO 2 (and hence temperature rise) is limited.
The differential sensitivity of surface ocean saturation state and pH to nCIE onset duration (figure S8) that occurs in our experiments has been noted previously (Hönisch et al 2012) and is a further example of how the timescale of carbon release interacts with the timescales of negative inorganic carbon cycle feedbacks (Hönisch et al 2012). A consequence of the weathering feedback is that when carbon input is sufficiently slow, the dominant sedimentary response is an increase in CaCO 3 accumulation rate. In our experiments, modeled sedimentary CaCO 3 concentrations drop initially in response to carbon addition, but after approximately 20-30 kyr (or less, for the shortest onset experiments) global sedimentary wt% CaCO 3 starts increasing. In fact, for most of our experiments this carbonate overshoot (e.g. Penman et al 2016) exceeds initial dissolution in terms of the change in global mean wt% CaCO 3 . It is important to note that we employ a fixed rain ratio (i.e. constant ratio of particulate organic carbon to particulate inorganic carbon throughout our experiments) and fixed remineralization profiles so that carbonate fluxes only change along with changes in calculated export productivity. The duration of carbonate dissolution is hence related to the weathering timescale and not to the duration of carbon input. As enhanced weathering delivers excess alkalinity to the oceans, carbonate burial will recover and rebound, regardless of whether carbon input to the atmosphere has ceased.
Comparison with past nCIEs A primary goal of our ensemble is to provide a framework for the interpretation of past geological events. Below, we provide examples in comparing our model experiments to three events: ETM-2, the PETM and OAE-1a.
ETM-2-The Eocene Thermal Maximum 2 (ETM-2, ∼54 Ma) is a relatively short-lived warming event marked by a symmetric nCIE of 1.0 to 1.5‰ and widespread deep sea CaCO 3 dissolution (Lourens et al 2005, Stap et al 2009 (figure 2(c)). The nCIE scenario modeled in this study that is most similar to the ETM-2 is the symmetric 1‰ nCIE with a total duration of 100 kyr without under-or overshoot in surface DIC δ 13 C ( figure 6(a)). To produce a nCIE of this size with an isotopically heavy carbon input (−6‰), approximately 5200 Pg C is required, resulting in atmospheric CO 2 increase of 900 ppm and SST increase of 2.6°C. An input of isotopically lighter carbon (−12 to −22‰) over 50 kyr produces a 1‰ nCIE with 2700 to 1500 Pg C respectively, and leads to a global mean SST increase of 1.5°C-0.8°C. Previous studies based on proxy data estimate (regional) SST increases during the ETM-2 of 3°C-5°C (Lourens et al 2005, which fits best to a carbon input of mantlederived volcanic source. However, the orbital pacing of a sequence of early Eocene hyperthermals, including the ETM-2, seems inconsistent with a volcanic origin. Increases in SST of 1.5°C-0.8°C produced by −12 and −22‰ carbon input are low compared to available SST proxies, but the maximum ETM-2 magnitude of 1.5‰ recorded in bulk carbonate is 0.5‰ larger than the 1‰ nCIE we compare to here and the true nCIE magnitude may have been even larger ). An nCIE of 1‰ driven by isotopically light (−60‰) methane requires ∼600 Pg C, produces an increase in atmospheric CO 2 of only 80 ppm and a 0.3°C SST increase. The latter would be nearly undetectable in most noisy proxy data. Furthermore, the widespread deep-sea carbonate dissolution associated with the ETM-2 does not result from experiments with isotopically lighter carbon sources. Less than 2.5 wt% global average CaCO 3 decline results from a 1‰ nCIE driven by either −60‰ or −22‰ carbon inputs. Hence, a larger carbon input (and therefore a carbon input with mean δ 13 C signature heavier than −22‰) is most consistent with the extent of environmental changes during the ETM-2. An alternative scenario that could explain the large temperature increase and globally widespread carbonate dissolution that is that the onset occurred on timescales much shorter than 50 kyr. nCIEs with an onset duration shorter than 50 kyr result in a significant drop in surface ocean calcite saturation state (Ω cal ) which promotes the dissolution of CaCO 3 (figures 6(b) and S8, S10). Likewise, the SST rise associated with a more rapid nCIE onset is larger compared to a nCIE onset duration over 50 kyr (figures 6(b) and S9).
PETM (∼56 Ma) is the greenhouse gas-driven global warming event most often suggested as an analog for the modern (e.g. Zeebe and Zachos 2013). Many more records of the PETM nCIE from a greater variety of settings are available in comparison to other Mesozoic and Cenozoic hyperthermals (McInerney and Wing 2011). The nCIE varies in size between different reservoirs, but the average bulk carbonate nCIE is closest in size to the 3‰ nCIEs modeled here (figures 2(b)). The PETM is well known for its asymmetric shape, with a rapid onset of <10 kyr and a total duration of ∼200 kyr (McInerney and Wing 2011, . None of our experiments shows an onset as rapid as the PETM-the most rapid modeled onset for a 100 kyr nCIE is 25 kyr. We compare the PETM to the 100 kyr, rapid onset 3‰ nCIE without under-or overshoot in surface DIC δ 13 C ( figure S11(b)). Estimated total carbon input is between 1600 to more than 40 000 Pg C depending on carbon source. Only the −12‰ carbon source leads to an increase in SST (5°C) consistent with compiled temperature proxy records (Dunkley Jones et al 2013), resulting in an increase in atmospheric CO 2 levels by ∼2500 ppm and declines in surface ocean pH and calcite saturation state (Ω cal ) of 0.47 and 2.7, respectively. The total carbon input for this scenario (>10 000 Pg C) is also consistent with recent estimates (Gutjahr et al 2017). However, the onset duration of the PETM is estimated to be shorter than 25 kyr and we can use our framework to explore the environmental consequences of more rapid onset duration as well. When a 3‰ nCIE onset is produced over 12.5 kyr (from the experiment with a modeled total duration of 50 kyr) using a carbon source with δ 13 C of −22‰, approximately 10 000 Pg carbon input is required. This increases atmospheric CO 2 levels by 3200 ppm, the global mean SST rises by 5.6°C, surface pH drops 0.57 units and Ω cal declines by 3.4.
OAE-1a-The Aptian Ocean Anoxic Event (OAE-1a) is one of the two major Cretaceous OAEs with a recorded nCIE preceding wide-spread organic carbon burial (Menegatti et al 1998, Herrle et al 2004, Jenkyns 2010) and might hold important information on carbon reservoir changes involved in the expansion of marine anoxia during this time. Aptian nCIEs are most commonly recorded with sizes ranging from −1 to −3‰. Using a nCIE at the upper limit of that interval from carbonate deposits at the Resolution Guyot (Jenkyns 1995) and age models based on cyclostratigraphy (Malinverno et al 2010), the nCIE at the onset of OAE 1a most closely resembles our 300 kyr-long simulations with rapid onset and δ 13 C values ending in an overshoot ( figure S11(a)). nCIEs of this shape and duration can be produced by biogenic methane emissions (δ 13 C−60‰) without drastically perturbing climate and ocean carbonate chemistry, but such a forcing mechanism cannot account for the reconstructed simultaneous rise in SST of 2°C-5°C (Mutterlose et al 2014, Naafs and, the atmospheric CO 2 increase of a few hundred to several thousand ppm , or the crisis of marine calcifying nannoplankton . A purely volcanic source of carbon results in a pCO 2 rise of 10-20 times pre-nCIE values, a stark decrease in surface pH and ocean calcite saturation, and a rise in global mean SST up to 10°C. Those environmental impacts are far more extreme than the proxy-based reconstructions suggest. Based on our framework and the age model provided by Malinverno (2010), the Aptian nCIE was most likely caused by a combination of organic and mantle carbon sources, resulting in the net emission of carbon of intermediate isotopic composition. However, recently published high-resolution records from shallow carbonate platforms suggest that the nCIE onset might have lasted significantly longer (Kuhnt et al 2011, Lorenzen et al 2013, Graziano and Raspini 2018, in which case a forcing from δ 13 C −6‰ volcanic carbon becomes more plausible for the OAE1a since longer nCIE onsets result in reduced impacts in climate and ocean chemistry, but a slight increase in the total mass of carbon that needs to be added to the ocean-atmosphere system (figure S10).

Summary/Conclusion
Our large ensemble of modeled nCIEs provides a template for interpreting past geologic events in terms of required carbon input and likely environmental consequences. We include the effects of various feedbacks such as changes in ocean solubility, carbon speciation, carbonate compensation and temperaturedependent weathering on the global carbon cycle. The earth system modeling approach allows us to track changes in the partitioning and fractionation of carbon between various carbon pools within the exchangeable reservoir. Our four dimensions of variability in defining modeled nCIEs (including size, carbon source, duration, and shape) adds two new dimensions compared to the conventional isotopic mass balance approach that neglects variability in duration and shape. Our results highlight the significant variability these two factors can generate in both the total carbon input required to generate a nCIE of a particular magnitude with a particular carbon source, but also the maximum rate of carbon input as a result of the interaction with modeled carbon cycle feedbacks. Thus, our experiments can be used to bracket the carbon emissions scenarios capable of driving any number of past nCIEs. Furthermore, these experiments demonstrate how environmental impacts scale with all four dimensions of variability. Hence, a comparison of available proxy data (e.g. changes in temperature, surface ocean pH, or changes in saturation state inferred from sedimentary carbonate dissolution) to our results can also be used to evaluate the likelihood of various carbon sources for a nCIE of a given magnitude.

Acknowledgments
While working on this manuscript PV and SKT were supported by a Heising-Simons Foundation award, MA was funded by NERC GW4+ PhD studentship S136361, and SEG was funded by NERC Independent Research Fellowship NE/L011050/1 and NERC large grant NE/P01903X/1. We thank James Rae and an anonymous reviewer for suggestions that greatly improved the manuscript.

Data availability
The specific version of the cGENIE.muffin model used in this paper is tagged as release v0.9.6 and has been assigned a DOI (https://doi.org/10.5281/ zenodo.3338584). The code is hosted on GitHub and can be obtained by cloning: https://github.com/ derpycode/cgenie.muffin, changing directory to cgenie.muffin, and then checking out the specific release: $ git checkout v0.9.6. Configuration files for the specific experiments presented in this paper can be found in the directory: genie-userconfigs\MS\vervoortetal.2019. Details of the experiments, plus the command line needed to run each one, are given in the README.txt file in that directory. All other configuration files and boundary conditions are provided as part of the release. A manual detailing code installation, basic model configuration, plus an extensive series of tutorials covering various aspects of muffin capability, experimental design, and results output and processing, is provided on GitHub. The latex source of the manual, along with a pre-built PDF file can be obtained, by cloning: https://github.com/derpycode/muffindoc. Model output files of this study are available from the corresponding author upon reasonable request.