the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 07 Aug 2026
Modeling the response of the marine carbon cycle to extreme CO2 injection events
Abstract. We explore how the response of a conceptual model of the marine carbon cycle depends on the way in which carbon is injected from the atmosphere. We find that, for single-injection pulses, the threshold amount required for a large response of the excitable system depends on pulse duration but not on its specific form. We do, however, see differences in the number of large transient responses in carbon and, correspondingly, the duration of the response for different pulse shapes. These differences are magnified as the system is pushed towards increased excitability and can be understood in terms of the geometry of an increasingly winding heteroclinic orbit. Inspired by Large Igneous Provinces (LIPs), we also consider random sequences of injection pulses. We find a wide range of possible responses for a given overall amount of injected carbon and duration, depending on the mean characteristics of the individual pulses. We also identify a resonance-like "Goldilocks" zone, in which intermediate pulse durations or arrival frequencies produce the largest number of repeated transients, and we test the framework with illustrative scenarios motivated by the Siberian Traps and Columbia River Basalt Group.
Status: open (extended)
-
RC1: 'Comment on egusphere-2026-4694', Anonymous Referee #1, 04 Sep 2026
reply
-
AC1: 'Reply on RC1', Punit Gandhi, 09 Oct 2026
reply
We thank the referee for the careful review and constructive comments that will help us improve the quality of our work. If given an opportunity to submit a revised version of the manuscript for ESD, we will take care to present our findings in a more consistent, self-contained way that focuses on the main novel results. In particular, we will be sure to standardize notation, provide more complete discussion of the model and remove introductory sections which may distract from the main focus of the paper.
In this manuscript we focus on exploring the extent to which details of CO2 injection from LIPs impact carbon cycle disruptions within the modeling framework of Rothman (2019). Because available data does not often tightly constrain the magnitude, duration and frequency of individual injection during these events of extreme volcanic activity, we aimed to explore the range of possible responses within the model for the same mean injection rate. We appreciate that the referee found the associated result with random injection pulses innovative. We recognize that the readability of the manuscript may have hindered clarity of our other novel results. Our revised manuscript will highlight our novel results more clearly, which are, briefly:
- Exploration of the impact that differing styles of CO2 injection can have on the response of the carbon system – i.e., not just the square pulse of Rothman, but explicitly investigating differences between square, ramp, exponential decay styles; these are all reasonable hypotheses for how LIPs may inject CO2 into Earth’s atmosphere.
- Exploration of the influence of random vs. steady pulses of CO2 injection; this expands the model to consider real world behavior of LIP eruptions.
- Exploration of pulse size and timing demonstrate that these exert a primary control on magnitude and duration of CO2 responses to LIP volcanism.
- Demonstration that slow nor fast, large nor small, high nor low frequency CO2 injection pulses alone produce the most significant climate perturbations. Instead, intermediate rate, magnitude, and pulse frequency have the largest climate response (anti-Goldilocks zone).
- Connection of modeled LIP injection and carbon cycle responses to real world data in the geologic record. This brings our study out of the theoretical and provides potential future directions of research to improve our understanding of the influence of different LIP eruptive styles on Earth’s climate and biosphere.
As pointed out by the referee, we do not undertake a rigorous dynamical systems analysis to establish the connection between the windings of the heteroclinic orbit and the observed response patterns. We are interested in carrying out such mathematical analyses of this phenomenon in future work, but believe that this falls outside of the scope of the current manuscript. For our exploration of the possible carbon cycle responses connected to LIPs, we consider the numerical simulations to provide sufficiently compelling evidence for the parameter regimes examined here. In a revised manuscript, we will be sure to distinguish our numerically-motivated interpretation from a rigorously demonstrated resonance mechanism.
We appreciate the referee’s feedback, and thank them for taking the time to provide us with specific comments and technical corrections that will improve the quality of our work and its presentation. We will be sure to incorporate these line-by-line comments and suggestions in a revised manuscript.
Citation: https://doi.org/10.5194/egusphere-2026-4694-AC1
-
AC1: 'Reply on RC1', Punit Gandhi, 09 Oct 2026
reply
-
RC2: 'Comment on egusphere-2026-4694', Anonymous Referee #2, 08 Sep 2026
reply
Publisher’s note: a supplement was added to this comment on 11 September 2026.
Gandhi et al use an existing conceptual model by Rothman to evaluate to what extent the rate and time of (LIP-like) C emissions affect carbon cycle feedback responses. The model consists of two feedbacks coupled to a simplified ocean with a (DIC) dampening timescale for air-sea exchange and ocean mixing. In certain parameter space, bistability exists which produces a non-proportional response to emissions when a threshold is met. Gandhi et al interrogate how the threshold changes with emissions rates with an ensemble of varying idealized and LIP-inspired emission scenarios.
General recommendation:
Rothman 2019 already demonstrated (by interrogating equations) that fast emission rates can lead to excitation of the system, amplifying the system response. Gandhi et al illustrate the implications with tested examples, i.e. a parameter sweep in mass-time space, this is hence a useful supplement to Rothman’s study. The methods and use of the model are sound, but the structure of the text could be improved. Authors mostly (not always) use appropriate citations. Figures are generally informative but somewhat difficult to interpret without context. The idealized scenarios to demonstrate patterns are useful but the LIP-inspired scenarios need improved justification. I believe the study fits the scope of ESD and contributes new insights though some critical points [see attached for full document] should be addressed, mainly improving readability and providing more, or better integrated, real-world application context.
-
AC2: 'Reply on RC2', Punit Gandhi, 09 Oct 2026
reply
We thank the referee for the constructive feedback, and very much appreciate the climate and biogeochemistry perspective. We find the comments and suggestions provided here to complement those of the other referee. In trying to straddle the line between an audience of applied mathematicians and Earth scientists, we missed the mark in some regards for both. If given the opportunity to submit a revision for ESD, we will incorporate the suggestions here to produce a more consistent, self-contained, and accessible text. In particular, we will be sure to include a more detailed introduction to the model, clearly define all terminology used, make figures easier to interpret, and try to make direct connections to Earth system processes throughout (as opposed to saving this to the end of the manuscript).
In the following, we reply to select significant questions raised by the referee to help clarify aspects of the work. We would be sure to address all more detailed suggestions and comments aimed at improving the presentation of our work in a revised manuscript.
Comment: 3. Testing of the square wave, exponential decay, and random perturbations is useful. Particularly the latter is interesting to explore how stochastic forcing affects the system and has direct links to LIPs. However, the sections need justifications for the choice of shape and magnitude of emissions, linking back to LIP-relevant processes. Needs citations and explanation about the type of LIP processes discusses here. E.g. extrusive, intrusive, thermogenic, explosive – all yield very different estimates for total mass-rate of injection. This is needed to support the choice of Mtot and distributions m̄, Tr. Also define ‘arrival time’. Same for Section 8, justify the choice for a 200,000 Gt of CO2 mass injection. This seems high – most papers suggest 30,000 Gt C (=110,000 Gt of CO2) for Siberian traps. The references cited all relate to the duration of emissions but not Mtot. Would be useful to explain what kind of volcanic processes are considered and their expected evolution over time to justify Tsqp and T.
Response: We agree that the geological motivation requires clearer justification. The pulse shapes are idealized forcing histories, rather than reconstructions of individual eruptions or intrusive episodes. In a revised manuscript we will distinguish geological constraints on carbon release and event duration from assumptions about pulse sizes, durations and arrival times. We will also verify the mass convention used in the calculations and report masses consistently. For reference, 200,000 Gt CO2 corresponds to approximately 54,500 Gt C, and this is consistent with, for example, table 1 of Jiang et al. (2022), PNAS 119(31). The source and time interval supporting the adopted Siberian Traps budget will be stated explicitly; a value used only to explore model sensitivity will be identified as such. We will distinguish total atmospheric emissions from effective input to the modeled ocean reservoir, and an ensemble budget fixed in expectation from one constrained in every realization.
Comment: 4. Related to point (1), describing the Earth system dynamics (how does the model respond in terms of the CaCO3 burial and surface respiration to LIP perturbation) is necessary to provide context and for palaeoclimatologists to determine whether these are reasonable responses within bounds of geological evidence. This explanation of model behaviour is lacking. Similarly, the test of what LIP-pulses excite the system in Section 8-9 describes the factual outcome of equations in terms of DIC peaks and recovery time but lacks the important discussion on paleoclimate implications. E.g. what does the mean number of DIC peaks mean? Is there geological evidence for these DIC peaks and strong variability?
Response: We agree that the response needs to be explained in terms of the modeled feedbacks. Carbon input raises DIC and reduces carbonate availability. In the model, increased shallow respiration amplifies the DIC disturbance, while reduced carbonate burial retains alkalinity and contributes to recovery, together with the prescribed relaxation term. In a revised manuscript, we will explain this sequence alongside the results and link this to broader climate implications. The mean peak count describes the number of repeated large transient responses across realizations; it is not a direct measure of extinction severity. We will discuss evidence for repeated carbon-cycle disturbances while distinguishing DIC peaks related to large transient responses in the model from carbon-isotope excursions and atmospheric CO2 injections, whose interpretation also depends on carbon sources, burial and mixing.
Comment: p. 9. Section 3.2. If I got this right, Ntr is the number of cycles before (c*,w*) is reached near equilibrium? Is Ntr equal to Nwind, does this relate to Fig. 3,? If so, reference. If not, illustrate elsewhere what is meant. (also clarify near p.16 Section 6).
Response: Ntr is the number of large DIC peaks counted as the system relaxes towards equilibrium after a perturbation, and this depends on the size of the perturbation. Nwind depends on system parameters, Figure 3 shows its dependence on Rtheta, and our numerical results for the isolated step perturbations considered here are consistent with Nwind acting as an upper bound on Ntr. No matter how large of a shift in nu (or equivalently, a downward shift in the initial condition of w), you can’t ever get more that Nwind peaks in the response. This comparison does not bound the cumulative number of peaks under repeated forcing. In a revised manuscript, we will clarify the scope and the criterion used to count a large peak.
Comment: p.10. Section 3.7 This seems to be a critical component of the manuscript but lacks paleoclimate background/implications. Can you provide some ranges of the timescales of ‘long forcing’ and ‘short pulses’? Are we talking yrs, kyr, Myrs? I realize the timescales depend on the processes and parameterization considered; identifying those timescales in the context of the two-feedback model use in this paper would be useful.
Response: In the baseline parameterization, the DIC relaxation time is 10,000 years. We use short and long forcing relative to this scale, rather than as absolute categories for LIP activity. In a revised manuscript, we will give dimensional durations alongside the forcing parameters and distinguish an individual pulse duration from the total duration of a pulse sequence.
Comment: Fig. 8. Interesting how a random distribution of carbon injections leads to seemingly periodic excitations (middle panels). Can you clarify whether the number of transient responses reflects the number of oscillations (e.g. as in Nwind), or if it reflects the number of times the injection pulse exceeds the threshold to excite the system, or a combination of both? Is there a consistent periodicity in this? If so, it could have interesting implications for the interpretation of periodicity in sedimentary sequences. The geologic significance of repeated excitations should be expressed somewhere.
Response: The number of transient responses is computed by counting large DIC peaks across an entire time series, which can include multiple CO2 injection pulses. This is different from (but related to) Nwind, which our numerical simulations suggest provides an upper bound on the number of transient responses to a single CO2 perturbation. A sufficiently strong CO2 injection pulse can initiate a sequence of up to Nwind large transient DIC responses, with the response depending on the state at which the pulse arrives. In the examples shown, additional injection pulses that arrive once a transient response is already underway have relatively little effect on the total number of transient responses experienced by the system. However, if there is another sufficiently strong injection pulse after the transient responses end, it can induce another sequence of transient responses. The near-regular spacing of the responses in the lower panel of figure 8B (Nwind=2 in this case) is consistent with the interpretation that above threshold injections are frequent enough that there isn’t much wait time after one sequence of responses ends and another one is initiated. We will make the distinction (and connection) between the number of large transient responses and Nwind more clear in a revised manuscript. We will also add a more detailed discussion of the near periodicity in response that can occur with random injections.
Comment: Fig. 9. I’m a little suspicious of DIC jumps that exceed 100,000 Gt. In reality, the DIC reservoir itself is unlikely reach such size as the carbonate ions needed to turn CO2 in HCO3 would run out. The max DIC jump may be mitigated when considering climate- dependent rock weathering (if that is not yet included...unclear). This is perhaps a limitation of the model that should be discussed.
Response: We agree that these extreme outputs need to be interpreted cautiously. They describe the reduced model under large forcing and should not be treated as quantitative reconstructions of past ocean DIC. Depletion of carbonate ions does not by itself impose a hard cap on total DIC, which also includes bicarbonate and dissolved CO2. Nevertheless, carbonate speciation, alkalinity balance, exchanges with other reservoirs and weathering feedbacks matter when assessing the physical plausibility of such large excursions. In a revised manuscript, we will clarify these limitations explicitly and distinguish the prescribed relaxation term from an explicitly modeled climate-dependent weathering feedback.
Comment: p.20. Note that at 250 Ma, the ‘ballasting hypothesis’ on which the R function is based (I believe?) may have operated differently due to limited pelagic calcifiers (Ridgwell, A., 2005. Marine Geology, 217(3-4)). Another way of changing R is via temperature (higher CO2=higher temp=more respiration). This essentially has a similar effect to the ballast hypothesis (higher CO2=less CaCO3 ballast =more respiration) so I don’t think it will affect your general conclusions, but the shape of the function likely differs depending on the assumed process with consequences for the threshold limits. This could be discussed or briefly mentioned for context.
Response: We agree that applying the ballast-based feedback to the end-Permian ocean requires qualification, given the different pelagic calcifier community. In a revised manuscript, we will clarify these differences with reference to Ridgwell (2005). Temperature-dependent respiration offers another possible amplification mechanism, but it is not represented explicitly in the present model. Similar feedback signs do not establish identical thresholds or dynamics, so we will avoid assuming that the quantitative results are unchanged under a different respiration law in our revision.
Comment: p.27 “Goldilocks zone” implies balance and a preferred or ideal state. Not sure this is the correct terminology here? Clarify what the Goldilocks zone means in this context.
Response: We will replace the “Goldilocks zone” expression with “anti-Goldilocks zone” to describe the intermediate range of pulse parameters that maximizes the mean number of large DIC excursions. We have seen this term used for analogous situations, in the context of economics and public policy, in order to highlight that the middle ground provides a worst-case combination of conditions.
Citation: https://doi.org/10.5194/egusphere-2026-4694-AC2
-
AC2: 'Reply on RC2', Punit Gandhi, 09 Oct 2026
reply
Viewed
Since the preprint corresponding to this journal article was posted outside of Copernicus Publications, the preprint-related metrics are limited to HTML views.
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 163 | 0 | 2 | 165 | 0 | 0 |
- HTML: 163
- PDF: 0
- XML: 2
- Total: 165
- BibTeX: 0
- EndNote: 0
Viewed (geographical distribution)
Since the preprint corresponding to this journal article was posted outside of Copernicus Publications, the preprint-related metrics are limited to HTML views.
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
General comments
The manuscript "Modeling the response of the marine carbon cycle to extreme CO2 injection events" by Gandhi et al. reviews a conceptual model of the marine carbon cycle originally formulated by Rothman (2017,2019), and studies the response of this excitable system to forcing pulses. The underlying motivation is extreme volcanic activity known as Large Igneous Provinces and its impact on the planet. The topic is scientifically relevant and within the scope of ESD. The authors review the model introduced by Rothman and discuss how three types of pulse perturbations affect the response of this excitable system. We note that the transient behaviour of the model against square injection pulses (one of the three types of perturbations considered in this work) had already been treated by Rothman (2019) as the authors correctly mention in their manuscript. The most (maybe only) innovative point seems to be the discussion of repeated random perturbations for a cumulative mass of injected CO2 showcasing three different dynamical behaviours depending on the mean size of the injection pulses, which also affects the injection frequency. The authors suggest to use the number of windings of the heteroclitic orbit as the measure of response of the system bringing into the discussion resonances and Goldilocks zones. The conclusions reached are of qualitative nature and do not go beyond the discussion of some numerical simulations. In other words, no rigorous dynamical analysis of the system is carried out beyond reviewing the autonomous case already presented by Rothman. The simulations are thoroughly described and therefore reproducible. Previous work is correctly referenced throughout. The title reflects the content of the paper and the abstract provides a complete and concise summary.
Specific comments
In my opinion the paper looks written a bit in a hurry, more similar to a collection of notes and results rather than trying to enhance the clarity for the reader's benefit. For example, I encourage the authors to (at least) review the following points:
1) the different notations for derivatives df/dt, f' and \dot f are interchangeably used, to the detriment of the manuscript's clarity.
2) The model is introduced rather quickly without providing enough context to understand its dynamics beyond referring to Rothman 2019 and plotting a bifurcation diagram. I encourage the authors to provide a more detailed presentation, possibly in appendix if they feel that this would disrupt the flow of the text. The same holds true for the reduction and its equilibrium (c*,w*), heteroclinic orbit, limit cycles, which are presented without any explanation.
3) On multiple occasions, the authors literally quote Rothman (2019), sometimes out of the blue, without any context (see for example third sentence after Figure 1, the end of page 9, the first sentence of Section 3.5, Section 3.7). The reader is inevitably forced to study Rothamn (2019) before being able of reading this manuscript.
4) Reading the introduction, I had the impression that the most innovative accomplishment of this manuscript was to discuss bifurcation-,noise-induced- and rate-induced tipping for the model by Rothman. However, Sections 3.5 and 3.6 seem completely superfluous as they are limited to few sentences directing the reader towards other papers and the phenomena of bifurcation- and noise-induced-tipping are never mentioned again in the text. Also Sections 3.7 and 4 on rate-induced tipping are respectively limited to reporting previous statements by Rothman, and setting up the ramp function and discussing one experiment in Figure 4.
Since seven people are listed as authors of this manuscript, I encourage a collective discussion on how to privilege the reader's perspective and make the paper clearer and more self-contained.
Technical corrections
- The constant C_p is used but not introduced in Section 2. The definition appears only in Table 1 without any discussion.
- At the end of page 3 the acronym DIC is used before being introduced in the following line
- At the end of page 4 the role of the variable c in determining an equilibrium is discussed before introducing the variable c.
- After formula (2b) it would be beneficial to write "where, for $\alpha>0$ we define," so that it is clear that in the ensuing formulas alpha is an exponent and not a super index. Otherwise there might be confusion on what S_\beta and \overline S_\gamma are.
-First line of subsection 3.3. The authors might want to write "In the (R_\theta,C_X) plane as this is consistent with the notation of Figure 2(b)
- End of page 9. The authors write "for the slice c_x=0.5 in panel a". However, it is not clear to which Figure they are referring to. If Figure 2, there is no c_x in panel a.
- Last sentence of Section 2.1. "summarized IN Table 1"