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: final response (author comments only)
-
RC1: 'Comment on egusphere-2026-4694', Anonymous Referee #1, 04 Sep 2026
-
AC1: 'Reply on RC1', Punit Gandhi, 09 Oct 2026
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
-
RC2: 'Comment on egusphere-2026-4694', Anonymous Referee #2, 08 Sep 2026
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
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
-
RC3: 'Comment on egusphere-2026-4694', Anonymous Referee #3, 12 Oct 2026
Using Rothman’s conceptual model of the marine carbon cycle, Gandhi et al. show that the timing of CO₂ release, not just the total amount, shapes how the ocean’s carbon system responds. The threshold for triggering large carbon excursions depends on pulse size and duration but not pulse shape. The number of repeated excursions, and so the length of the disruption, depends sensitively on timing: an abrupt pulse that ends near the peak of an excursion triggers additional ones, and under random pulse sequences the count peaks at intermediate pulse sizes and frequencies. Notably, the duration of disruption and its peak severity are controlled by different features of the forcing. Applied to the Siberian Traps and Columbia River Basalts, the model suggests that a more easily disturbed greenhouse ocean could explain why the Siberian Traps were far deadlier, though that comparison rests on loosely constrained assumptions.
The model itself isn’t new. All equations and parameters come from Rothman (2017, 2019), and the threshold behavior for square pulses was already worked out there. What is new is the systematic exploration of forcing structure. The paper compares pulse shapes, maps the number of excursions rather than just whether one occurs, links that count to the geometry of the model’s underlying dynamics, and introduces random pulse sequences inspired by how LIPs actually erupt in episodes. In short, this is an incremental but meaningful extension of an existing model, whose main contribution is a new way of thinking about eruption timing rather than new physics or new data.
Major Comments.
- The manuscript as a whole felt choppy. I think this could be fixed with a restructuring or combining of subsections and more detail explaining the relationships or purpose from one idea to the next.
- From a dynamical-systems standpoint, essentially everything structural comes from Rothman: the excitable regime, the bifurcations, the reduced saddle and heteroclinic connection, and the rate–mass dichotomy. What this paper adds is numerical simulation under new forcing protocols, with the results interpreted through that existing geometry. There are no new theorems, no new bifurcation analysis, and no reduction that explains the observed structure. The winding count in Figure 3 is computed numerically, the resonance mechanism is described qualitatively from phase portraits, and the random-forcing results are Monte Carlo averages with no stochastic analysis behind them. I think I would make it clearer to the reader from the start that the novelty of this manuscript is that it runs new forcing experiments on a known system, with the geological framing as the main motivation, rather than contributing anything new to dynamical systems itself. I think (potentially) why there is the disconnect is that the paper's introduction talks about bifurcation-, rate-, and noise-induced tipping, as well as the first section recaps Rothman's model and dynamical systems derivation, and I thought we were going to build on that.
- "For each case, 100 trials are simulated...": how did you decide you had enough simulations to get an accurate picture? Did you run Monte Carlo simulations until convergence? I think there needs to be more discussion on this.
- Several figures rely on color to differentiate parts of the plot and become illegible when printed in grayscale. Please revise figures to ensure they are adequately legible in grayscale. It definitely hindered interpretation. I also felt that the figures needed more explanation as I was not always sure what the aim was/there was a lot of information to parse out.
- I recognize that you need a solid understanding of Rothman's model and the dynamical system before building on it, but it felt like you spent a lot of time just writing what is already out there, but simultaneously you rushed the description/derivation as a first time reader would not necessarily know all of the jargon. I don't have a good suggestion for fixing this at the current time.
Technical comments.
- pg.1 giga-tons should be gigatons
- pg. 5, last sentence of section 2.1: “the conversions summarized Table 1” is missing “in.”
- pg.14 , Fig. 5 caption: it refers to a “white line in panel (a),” but I cannot see a white line.
- pg. 20, Section 8 intro: “the system 2” should be “system (2).”
- pg.25, Fig. 12: “Imnana” should be “Imnaha,” as spelled in Table 2.
- p. 29: “total CO₂ injection was estimated injection was estimated” has a duplicated phrase.
- p. 31: “complimentary pathways” should be “complementary.”
- In general, you change derivative notation without defining the definition of the dot notation. Further, you should be careful since you have not defined it as you rescale time and then use the same notation. Sloppy from the dynamical systems perspective.
- You used a lot of dynamical systems Jargon used without formal introduction which could be an issue for non DS readers. However, I think the most important would include:
- Heteroclinic orbit, which the whole paper relies on but never defines
- Hopf bifurcation
- Saddle-node of periodic orbits (SNPO). It’s named, but what it is isn’t explained
- First Lyapunov coefficient
- “Ghost” of a limit cycle
- Quasipotential, Kramers-type law, large-deviation, and first-exit/first-passage time
- It may be silly, but a quick def of fixed points being where deriv=0 wouldn't hurt.
Citation: https://doi.org/10.5194/egusphere-2026-4694-RC3
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 | |
|---|---|---|---|---|---|
| 206 | 0 | 2 | 208 | 0 | 0 |
- HTML: 206
- PDF: 0
- XML: 2
- Total: 208
- 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"