# Characteristic disruptions of an excitable carbon cycle

Edited by Donald E. Canfield, Institute of Biology and Nordic Center for Earth Evolution, University of Southern Denmark, Odense M, Denmark, and approved June 4, 2019 (received for review March 26, 2019)

## Significance

The great environmental disruptions of the geologic past remain enigmatic. Each one results in a temporary change in the oceans’ store of carbon. Although the causes remain controversial, these changes are typically interpreted as a proportionate response to an external input of carbon. This paper suggests instead that the magnitude of many disruptions is determined not by the strength of external stressors but rather by the carbon cycle’s intrinsic dynamics. Theory and observations indicate that characteristic disruptions are excited by carbon fluxes into the oceans that exceed a threshold. Similar excitations follow influxes that are either intense and brief or weak and long-lived, as long as they exceed the threshold. Mass extinction events are associated with influxes well above the threshold.

## Abstract

The history of the carbon cycle is punctuated by enigmatic transient changes in the ocean’s store of carbon. Mass extinction is always accompanied by such a disruption, but most disruptions are relatively benign. The less calamitous group exhibits a characteristic rate of change whereas greater surges accompany mass extinctions. To better understand these observations, I formulate and analyze a mathematical model that suggests that disruptions are initiated by perturbation of a permanently stable steady state beyond a threshold. The ensuing excitation exhibits the characteristic surge of real disruptions. In this view, the magnitude and timescale of the disruption are properties of the carbon cycle itself rather than its perturbation. Surges associated with mass extinction, however, require additional inputs from external sources such as massive volcanism. Surges are excited when ${\text{CO}}_{2}$ enters the oceans at a flux that exceeds a threshold. The threshold depends on the duration of the injection. For injections lasting a time ${t}_{i}\gtrsim \mathrm{10,000}$ y in the modern carbon cycle, the threshold flux is constant; for smaller ${t}_{i}$, the threshold scales like ${t}_{i}^{-1}$. Consequently the unusually strong but geologically brief duration of modern anthropogenic oceanic ${\text{CO}}_{2}$ uptake is roughly equivalent, in terms of its potential to excite a major disruption, to relatively weak but longer-lived perturbations associated with massive volcanism in the geologic past.

### Sign up for PNAS alerts.

Get alerts for new articles, or get an alert when an article is cited.

## Introduction

Earth’s carbon cycle is a loop between photosynthesis, which converts carbon dioxide (${\text{CO}}_{2}$) to organic carbon, and respiration, which converts organic carbon back to ${\text{CO}}_{2}$ (1–4). The cycle has undergone many disruptions throughout Earth’s history (5–7). These events are expressed in the geologic record as relatively abrupt and large changes in the isotopic composition of sedimentary carbon (8) compared with background fluctuations. Fig. 1 shows 2 examples in one 800-ky time series. Such events are typically attributed to changes in the fluxes and concentrations of carbon (8–10). Thus, Fig. 1 could reflect the transient addition of 2 pulses of isotopically light carbon originating, for example, from the respiration of a previously unreactive reservoir of organic carbon (10, 11). Such disruptions have also been interpreted as geochemical responses to other sources of environmental change, including variations in Earth’s orbital parameters (12); dissociation of methane hydrate (13); bolide impacts (14); biogeochemical innovations (15, 16); and changes in chemical weathering (17), organic carbon burial (18), and volcanic emissions (19). These interpretations typically treat the marine carbon cycle as a passive recorder of an externally imposed stress. The response imprinted in the geochemical record is then implicitly assumed to be proportional to the magnitude of the forcing (8, 16, 20).

Fig. 1.

Given the variety of possible forcing mechanisms, it comes as no surprise that the size and timescale of disruptions vary immensely, by 2 orders of magnitude during the Phanerozoic (0 to 542 Ma) (7). Yet the same data show that major fluctuations of the carbon cycle exhibit a characteristic rate of change. Events associated with mass extinction tend to exceed the characteristic rate, whereas others appear bounded by it (7). It seems unlikely that a rich diversity of stressors, expressed as proportionate responses in the geochemical record, would exhibit such uniformity.

Here I formulate and analyze an elementary mathematical model that shows how the characteristic rate can instead emerge within the carbon cycle. The model portrays the upper ocean as a chemical reactor open to an incoming flux of dissolved calcium carbonate from rivers and an outgoing flux representing carbonate burial in sediments. Within the reactor, the concentrations of carbonate species respond not only to imbalances in inputs and outputs, but also to imbalances in the biological consumption and production of ${\text{CO}}_{2}$. This framework reveals a mechanism for autocatalytic amplification of a small but finite perturbation of a stable steady state followed by relaxation back to the same steady state. The process is analogous to the excitation of an action potential (nerve impulse) in a neuron (21). Here, excitations manifest themselves as transient ocean acidification; i.e., a temporary increase in the concentration of carbon dioxide in the upper ocean. This change is distinguished by a characteristic rate.

These results suggest that the magnitude of a carbon cycle’s disruption is determined not by the strength of the cycle’s perturbation but rather by the intrinsic dynamics of the system itself. Once the addition of ${\text{CO}}_{2}$ to the oceans passes a threshold, the rate of amplification and the eventual severity of the resulting environmental change should be independent of the detailed history of the perturbation. Moreover, the model indicates that the consequences of fast forcing at human timescales may be similar to the outcome of slow forcing at geologic timescales (7). Within this framework, the greater rates of change associated with mass extinction events are straightforwardly interpreted as the consequence of forcing sustained beyond the threshold, due, for example, to massive volcanism (22–25).

This paper is organized as follows. The first section presents a 2-component dynamical system that represents important features of the marine carbonate system. Next, this paper analyzes the stability of the model’s steady state. Three dynamical regimes receive focus: The limit cycle that results when the steady state is unstable, the excitations that result when the steady state is stable, and the model’s behavior when both the steady state and the limit cycle are stable. Particular attention is given to the case of excitations and their size and timescale. This paper then discusses the implications of these findings for interpreting the geochemical record of past disruptions and predicting the carbon cycle’s response to modern anthropogenic perturbations. The conclusion highlights strengths and weaknesses of this paper’s principal contributions and points the way toward future progress.

## Two-Component Dynamical System

This section constructs a mathematical model to help explore how instabilities might occur in the marine carbon cycle. Among the many potentially relevant feedback mechanisms (27), I focus on the interactions of shallow-ocean respiration with fluxes of carbon into and out of the shallow ocean. These feedbacks are expressed below as a dynamical system composed of 2 ordinary differential equations. The objective is to provide a rational context within which one can show how mechanisms operating within the carbon cycle can give rise to general features of its past disruptions. Such understanding follows from the identification of the ways in which the qualitative dynamics of the system—e.g., its stability and bifurcations—can derive from its intrinsic interactions. The theory of dynamical systems makes such connections possible (28–32).

Although these objectives require the identification of relevant mechanisms, they do not necessitate detailed parameterizations. The developments below nevertheless give significant attention to the way in which specific processes operate in the real carbon cycle. I begin by stating essential aspects of carbonate chemistry. The model is then stated in its most general form, after which pertinent mechanisms are expressed mathematically. A final change of variables then leads to a model suitable for classical stability analysis.

### The Carbonate System.

The model expresses the evolution of total alkalinity and the total dissolved inorganic carbon. Total alkalinity is defined by the sum (2, 33)Total alkalinity represents the excess of bases over acids. Its dominant contributions come from bicarbonate ions (${\text{HCO}}_{3}^{-}$) and carbonate ions (${\text{CO}}_{3}^{2-}$); the factor of 2 attached to the ${\text{CO}}_{3}^{2-}$ concentration derives from its double-negative charge (33). Total dissolved inorganic carbon (DIC) is symbolized by $w$, mnemonically recalled as the “whole” of inorganic carbon:The carbonate system is taken to always be in thermal equilibrium according to the reactions (2, 33)The 2 conserved quantities ($a$ and $w$) and the 2 reactions of Eq.

$$a=\left[{\text{HCO}}_{3}^{-}\right]+2\left[{\text{CO}}_{3}^{2-}\right]+\left[{\text{B}\left(\text{OH}\right)}_{4}^{-}\right]+\left[{\text{OH}}^{-}\right]-\left[{\text{H}}^{+}\right].$$

[1]

$$w=\left[{\text{CO}}_{2}\right]+\left[{\text{HCO}}_{3}^{-}\right]+\left[{\text{CO}}_{3}^{2-}\right].$$

[2]

$${\text{CO}}_{2}+{\text{H}}_{2}\text{O}\rightleftharpoons {\text{HCO}}_{3}^{-}+{\text{H}}^{+}\rightleftharpoons {\text{CO}}_{3}^{2-}+2{\text{H}}^{+}.$$

[3]

**3**specify the equilibrium of the carbonate system. There are 6 unknowns: Not only $a$ and $w$, but also the concentrations of ${\text{CO}}_{2}$, ${\text{CO}}_{3}^{2-}$, ${\text{HCO}}_{3}^{-}$, and ${\text{H}}^{+}$. Specification of any 2 of these quantities therefore provides the other 4 via equations (1–3, 33).### General Formulation.

I consider a well-mixed open system in which dissolved calcium carbonate (${\text{CaCO}}_{3}$) is delivered to the shallow ocean by rivers and output from the shallow ocean ultimately into sediments. The assumption of a well-mixed system implies that the model addresses timescales longer than the timescale of global ocean circulation [about 1,000 y (2, 3)]. The precise thickness of the shallow layer is unimportant; however, it must be small compared with the average depth of the oceans. This assumption allows exchanges of carbon between the shallow ocean and the deep ocean and sediments to be expressed in analogy to the exchange of heat between a small system and a large heat bath. In other words, the concentration of carbon fluctuates within the shallow layer but the mass of carbon below it is effectively constant. Thus, the model explicitly considers only the chemical state of the upper ocean.

Fig. 2 illustrates the model. Its mathematical form readswhere dots represent time derivatives. The concentrations $a$ and $w$ are expressed in units of $\mathrm{\mu}\text{mol}\cdot {\text{kg}}^{-1}$. The factor of 2 arises because 2 molar equivalents of alkalinity are associated with each mole of ${\text{CaCO}}_{3}$ (33). Expressing the model in terms of the conserved quantities $a$ and $w$ rather than, e.g., the concentrations of ${\text{CO}}_{2}$ and ${\text{CO}}_{3}^{2-}$ makes it possible to temporarily forestall a discussion of the changing chemical equilibrium.

$$\u0227=2\left[\hspace{0.17em}{j}_{\mathrm{i}\mathrm{n}}-B\left(a,w\right)\right]$$

[4]

$$\stackrel{\u0307}{w}=\left(1+\nu \right){j}_{\mathrm{i}\mathrm{n}}-B\left(a,w\right)+R\left(a,w\right)-\left(w-{w}_{0}\right)/{\tau}_{w},$$

[5]

Fig. 2.

The flux ${j}_{\mathrm{\text{in}}}$ represents the rate at which dissolved ${\text{CaCO}}_{3}$ is added to the oceans by rivers in terms of an equivalent change in the whole-ocean concentration of dissolved ${\text{CaCO}}_{3}$, expressed in units of $\mathrm{\mu}\text{mol}\cdot {\text{kg}}^{-1}\cdot {\text{y}}^{-1}$. The flux $B\ge 0$ denotes the rate at which ${\text{CaCO}}_{3}$ is output from the system, which, in the real ocean, corresponds to burial in sediments. However, the model makes no dynamical distinction between the deep sea and sediments. Thus, $B$ represents the rate at which ${\text{CaCO}}_{3}$ exits the shallow layer without returning. The dimensionless term $\nu $ determines the strength of an external source of ${\text{CO}}_{2}$, such as volcanic emissions.

The flux $R\ge 0$ represents the extent to which upper-ocean production of ${\text{CO}}_{2}$ by respiration exceeds a theoretical baseline. The small change in alkalinity associated with respiration (33) is neglected. When $\nu =0$, the baseline case $R=0$ is associated with a steady state in which $w={w}_{0}$ and the riverine influx of dissolved ${\text{CaCO}}_{3}$ is precisely balanced by burial (i.e., $B={j}_{\mathrm{\text{in}}}$). In this case, the model’s upper ocean would act merely as a way station for the temporary storage of dissolved carbonate. Positive $R$, on the other hand, acts to increase $w$ above ${w}_{0}$. The model assumes that this increase is resisted linearly, so that if the baseline $R=0$ were restored, $w$ would return to ${w}_{0}$ at the timescale ${\tau}_{w}$. At long timescales ($\gg 1{0}^{2}$ y) in the real ocean, the processes responsible for such dissipation include the riverine influx of dissolved ${\text{CaCO}}_{3}$ and the dissolution and precipitation of ${\text{CaCO}}_{3}$ on the seafloor (34, 35). These processes act as the ocean’s “homeostat” or “pHstat” at a wide range of timescales (2, 35). Aspects of the homeostat that act simultaneously on DIC and alkalinity respond to the difference between ${j}_{\mathrm{\text{in}}}$ and $B$. However, the homeostat must also influence DIC independently by opposing the respiration feedback; otherwise DIC would grow without bound when $R>0$. Eq.

**5**represents this negative feedback by the term proportional to $w-{w}_{0}$. The homeostat’s dominant characteristic timescale, denoted here by ${\tau}_{w}$, is about 10 ky (36–38) in the modern ocean.Eqs.

**4**and**5**describe the general form of a dynamical system in which the riverine influx of alkalinity and DIC produce internal feedbacks that collectively act to determine the extent to which alkalinity and DIC accumulate before exiting the system. The description of the model is completed by the functional specification of the feedbacks $B$ and $R$, to which we now turn.### The Burial and Respiration Fluxes.

The export and burial of ${\text{CaCO}}_{3}$ in the real ocean depends on the physical environment (e.g., temperature and pressure), the chemical environment, and ecological factors affecting biogenic carbonate precipitation (e.g., phytoplankton community composition). When seawater is undersaturated with carbonate ions, ${\text{CaCO}}_{3}$ dissolves; when it is supersaturated, ${\text{CaCO}}_{3}$ can precipitate (33). At low saturations, physical and biogenic precipitation—and therefore carbonate output—is impeded. The dependence of precipitation on pressure manifests itself as a critical depth beneath which carbonate dissolves. Based in part on earlier work (39), I express these dependencies bywhere the carbonate ion concentration $c\left(a,w\right)$ derives from the equilibrium of the carbonate system (33) andis the sigmoidal function shown in

$$B\left(a,w\right)=B\left[c\left(a,w\right)\right]=b{j}_{\mathrm{\text{in}}}s\left(c,{c}_{p}\right),$$

[6]

$$s\left(c,{c}_{p}\right)=\frac{{c}^{\gamma}}{{c}^{\gamma}+{c}_{p}^{\gamma}}$$

[7]

*SI Appendix*, Fig. S1. The exponent $\gamma $ parameterizes the sharpness of the transition, ${c}_{p}$ is the transitional carbonate ion concentration where $s=1/2$, and $b{j}_{\mathrm{\text{in}}}$ denotes the maximum rate of carbonate burial. The apparent requirement that $c$ be computed from $a$ and $w$ would appear to be a serious complication, but a change of variables discussed below solves this problem. The introduction of the exponent $\gamma $ may also appear unwelcome; however, this paper’s conclusions require only that the function $s$ be sigmoidal.*SI Appendix*, Fig. S1 suggests that $\gamma \simeq 4$ corresponds to the modern ocean.Respiration is also sensitive to environmental conditions. For example, warmer temperatures should result in increased respiration rates relative to production rates (27, 40–42). Possible consequences include an upward shift of respiration and a positive feedback for shallow-ocean and atmospheric ${\text{CO}}_{2}$ levels (27, 41, 43, 44). Alternatively, carbonate production by pelagic calcifiers may decrease at lower pH (45–48). The “ballast hypothesis” (49–51) suggests that the export of organic carbon out of the photic zone is facilitated by its association with biogenic carbonate. A decreased supply of ballast would decrease export of organic carbon and therefore increase respiration rates in the shallow ocean. A positive feedback is again possible (27, 52, 53).

Neither the ballast (54–58) nor the temperature (42, 59–61) feedback is free of ambiguity. However, both are reasonable hypotheses. This paper addresses the case of the ballast feedback, with the expectation that the lessons learned will inform consideration of other possible respiration feedbacks, including those involving temperature.

The essential idea is that lower concentrations of ${\text{CO}}_{3}^{2-}$ lead to less biogenic carbonate production and therefore less ballast, less associated organic carbon export, and more upper-ocean respiration (27, 45–53). These observations suggest that the respiration flux $R$ decreases with increasing carbonate ion concentration $c$. For $c$ below a low concentration of carbonate ions, $R$ should be near its maximum. Conversely, for $c$ above a high concentration, respiration in the shallow layer is assumed to reach its baseline $R=0$. The baseline case therefore corresponds to the upper-ocean respiration flux that would continue to occur despite elevated production of carbonate ballast.

This reasoning suggests that $R$ decreases with $c$ similar to the way $B$ increases with $c$ (where $\overline{s}=1-s$, $\theta {j}_{\mathrm{\text{in}}}$ is the maximum value of $R$, and ${c}_{x}$ is the “crossover” concentration that marks the halfway point in the transition from high to low $R$. As for the burial function of Eq.

*SI Appendix*, Fig. S1), but with a possibly different amplitude and transitional concentration. I therefore write$$R\left(a,w\right)=R\left[c\left(a,w\right)\right]=\theta {j}_{\mathrm{\text{in}}}\overline{s}\left(c,{c}_{x}\right),$$

[8]

**6**, the sharpness of the transition parameterized by $\gamma $ is a minor detail.*SI Appendix*, Fig. S2 illustrates how the burial and respiration fluxes of Eqs.**6**and**8**modify Fig. 2.### Transformation from $\u0227$ to $\u010b$.

*SI Appendix*shows that the rate of change of the ${\text{CO}}_{3}^{2-}$ concentration can be approximated by

$$\u010b=f\left(c\right)\left(\u0227-\stackrel{\u0307}{w}\right).$$

[9]

$$f\left(c\right)={f}_{0}\frac{{c}^{\beta}}{{c}^{\beta}+{c}_{f}^{\beta}}$$

[10]

*SI Appendix*, Fig. S3). Substitution of Eqs.

**6**and

**8**into Eqs.

**4**and

**5**and insertion of the resulting expressions into Eq.

**9**then remove the implicit dependence of $c$ on $a$ and $w$. The new system reads

$$\u010b/f\left(c\right)=\mu \left[1-bs\left(c,{c}_{p}\right)-\theta \overline{s}\left(c,{c}_{x}\right)-\nu \right]+w-{w}_{0}$$

[11]

$$\stackrel{\u0307}{w}=\mu \left[1-bs\left(c,{c}_{p}\right)+\theta \overline{s}\left(c,{c}_{x}\right)+\nu \right]-w+{w}_{0},$$

[12]

$$\mu ={j}_{\mathrm{\text{in}}}{\tau}_{w}$$

[13]

**11**and

**12**is henceforth called the carbon-cycle model. Its parameters are listed in

*SI Appendix*, Table S1. Of these, only $\mu $, $b$, $\theta $, ${c}_{x}$, ${c}_{p}$, and $\nu $ are of interest. The remaining 5 parameters are set to correspond to the equilibrium chemistry or the properties of the modern ocean. Although they are retained in the analysis that follows, they serve only to maintain realism.

### Steady State and Its Stability.

Steady states, or fixed points, occur where $\u010b=\stackrel{\u0307}{w}=0$. Under this condition, the addition of Eqs. which exists only for $b>1$. Substitution of ${c}^{*}$ for $c$ in Eq. Note that the untransformed model (Eqs.

**11**and**12**yields $bs\left(c,{c}_{p}\right)=1$. Substituting Eq.**7**and solving for $c$, one then finds the unique steady-state carbonate ion concentration$${c}^{*}={\left(b-1\right)}^{-1/\gamma}{c}_{p},$$

[14]

**12**then provides the steady-state DIC concentration$${w}^{*}={w}_{0}+\mu \left(\theta +\nu -\frac{\theta {c}_{p}^{\gamma}}{\left(b-1\right){c}_{x}^{\gamma}+{c}_{p}^{\gamma}}\right).$$

[15]

**4**and**5**) yields the same steady state after substitution of Eqs.**6**and**8**.*SI Appendix*shows that parameter values consistent with the modern ocean (*SI Appendix*, Table S1) predict values of ${c}^{*}$ and ${w}^{*}$ that are also consistent with the modern ocean.*SI Appendix*also derives and interprets the results of a linear stability analysis. The requirements for stability are expressed in terms of the trace $\text{Tr}$ and determinant $\mathrm{\Delta}$ of the Jacobian matrix of Eqs.

**11**and

**12**evaluated at the fixed point. These calculations yield

$$\text{Tr}=-1-\frac{\mathrm{\Delta}}{2}\left(1-\frac{\theta b{c}_{p}^{\gamma}{c}_{x}^{\gamma}}{{\left[\left(b-1\right){c}_{x}^{\gamma}+{c}_{p}^{\gamma}\right]}^{2}}\right),$$

[16]

$$\mathrm{\Delta}=\frac{2\mu \gamma {f}_{0}{\left(b-1\right)}^{\frac{1}{\gamma}+1}{c}_{p}^{\beta -1}}{b{\left(b-1\right)}^{\beta /\gamma}{c}_{f}^{\beta}+b{c}_{p}^{\beta}}.$$

[17]

Fig. 3.

## Limit Cycle

*SI Appendix*shows that the fixed point becomes unstable via a Hopf bifurcation, which leads to periodic nonlinear oscillations called limit cycles. Fig. 4 depicts a stable limit cycle in the phase plane of $c$ and $w$.

*SI Appendix*, Fig. S4 illustrates the time series $c\left(t\right)$ and $w\left(t\right)$ on the limit cycle, along with the time evolution of quantities derived from these 2 quantities.

Fig. 4.

To understand how the limit cycle operates, consider the evolution of the periodic trajectory from the point where $w$ is at its minimum, i.e., at the lower crossing of the nullcline $\stackrel{\u0307}{w}=0$. Here, the ballast feedback begins to increase DIC by the addition of ${\text{CO}}_{2}$, which acts to consume ${\text{CO}}_{3}^{2-}$ ions via the reaction of

*SI Appendix*, Eq.**S8**. However, the relative rate at which the ${\text{CO}}_{3}^{2-}$ concentration decreases compared with the rate at which DIC increases becomes progressively slower as the ocean becomes less buffered at lower ${\text{CO}}_{3}^{2-}$ concentrations. Eventually the ${\text{CO}}_{3}^{2-}$ concentration reaches a minimum where the loss of ${\text{CO}}_{3}^{2-}$ via burial and the ${\text{CO}}_{2}$ feedback is balanced by the homeostatic response, and it crosses the nullcline $\u010b=0$. However, DIC continues to rise until it is overcome by the homeostatic feedback, reaching a maximum where it crosses the nullcline $\stackrel{\u0307}{w}=0$. The homeostatic feedback then acts to decrease $w$ while increasing $c$. The ${\text{CO}}_{3}^{2-}$ concentration reaches its maximum when carbonate burial becomes stronger than the addition of dissolved ${\text{CaCO}}_{3}$ from rivers and the weakened homeostatic feedback. A rapid upward jump of the burial rate then acts to reduce both the ${\text{CO}}_{3}^{2-}$ and DIC concentrations. Eventually DIC reaches its minimum and the cycle repeats.A single period of the limit cycle embodies several attributes of the marine carbonate system that have attracted considerable attention in recent years. The rise of DIC, for example, is initially accompanied by a rise of ${\text{CO}}_{2}$ (

*SI Appendix*, Fig. S4*C*) and a reduction of pH (*SI Appendix*, Fig. S4*E*). This process is called ocean acidification (62). The acid is eventually neutralized by the dissolution of seafloor carbonate and other homeostatic processes, but only slowly, at a modern timescale of about 10 ky—part of the long legacy of any rise in ${\text{CO}}_{2}$ levels (34). Finally, the surfeit of dissolved ${\text{CaCO}}_{3}$ must eventually be buried. In the limit cycle pictured here, carbonate burial occurs as a sharp pulse (*SI Appendix*, Fig. S4*H*) analogous to the “cap-carbonate” deposits that follow global glaciations (63).Periodic disruptions of the carbon cycle independent of periodic forcing may have occurred in the geologic past (64–67). The following section shows, however, that the potential existence of a limit cycle is more important than its actual existence. Three aspects of the limit cycle are central to that discussion: Its amplitude, its period, and its dynamical origin.

Unlike linear oscillations, the amplitude or size of a limit cycle is independent of the initial condition. As shown in Fig. 4, trajectories originating inside the limit cycle grow outward, whereas trajectories originating outside move inward. The resulting amplitude is therefore a property of the system itself rather than its perturbation.

However, the size of the limit cycle does depend on the system’s parameters. An appropriate measure of size along the DIC axis is the difference $\mathrm{\Delta}w={w}_{\mathrm{\text{max}}}-{w}^{*}$, where ${w}_{\mathrm{\text{max}}}$ is the maximum value of $w$ in the cycle; one similarly has $\mathrm{\Delta}c={c}_{\mathrm{\text{max}}}-{c}^{*}$. where $A={\left(\mathrm{\Delta}w\mathrm{\Delta}c\right)}^{1/2}$ is a measure of the amplitude of the limit cycle in units of $\mathrm{\mu}\text{mol}\cdot {\text{kg}}^{-1}$ and ${T}_{0}\sim 1$ is the dimensionless period of an infinitesimal limit cycle. The factor of 2 reflects an assumption that roughly twice as much time is required to go through a complete cycle compared with the time it takes to change $w$ or $c$ by its typical excursion.

*SI Appendix*shows that $\mathrm{\Delta}w$ and $\mathrm{\Delta}c$ generally increase with $\mu $ and $\theta $ but decrease with $b$. The limit cycle’s size grows roughly linearly with $\mu $ because $\mu $ sets the rate of all processes other than homeostatic damping. Because rates and size scale similarly with $\mu $, the dimensionless period $T$ of all but the smallest limit cycles should be effectively independent of $\mu $. These observations imply that$$T\simeq {T}_{0}+2A/\mu ,$$

[18]

*SI Appendix*, Fig. S5 shows that Eq.

**18**is a reasonable approximation. It characterizes the limit cycle well partly because the complicated dependence of $A$ and $T$ on the system’s parameters other than $\mu $ (

*SI Appendix*) is similarly expressed by $A$ and $T$. The physical interpretation of Eq.

**18**is readily exposed by rewriting it in terms of dimensional time. Setting $\tau =T{\tau}_{w}/2$, ${\tau}_{0}={T}_{0}{\tau}_{w}/2$, and recalling $\mu ={j}_{\mathrm{\text{in}}}{\tau}_{w}$, one finds

$$A\simeq {j}_{\mathrm{\text{in}}}\left(\tau -{\tau}_{0}\right).$$

[19]

The manner in which the limit cycle arises is also of interest. There are 2 kinds of Hopf bifurcations—subcritical and supercritical (21, 28, 29).

*SI Appendix*performs a standard calculation (21, 28, 29) to determine where each Hopf bifurcation occurs here; the results are color coded along the curves of Fig. 3. In the carbon-cycle model, the subcritical Hopf bifurcation is adjacent to a bistable region where the stable fixed point and a stable limit cycle appear with an unstable limit cycle. The following section establishes the importance of this region for the interpretation of the model’s response to perturbations of the stable fixed point.## Excitability

Now suppose the system is in a stable steady state with no injection of ${\text{CO}}_{2}$ ($\nu =0$). At time $t=0$ we set $\nu =0.35$, which corresponds to injecting additional DIC into the oceans at a rate equal to 35% of riverine inputs. From Eqs.

**14**and**15**, we see that ${\text{CO}}_{2}$ injection increases ${w}^{*}$ by $\nu \mu $, but ${c}^{*}$ is unchanged. The new fixed point retains its stability, because stability is independent of $\nu $ (Eqs.**16**and**17**). However, the state of the system must now relax to the new fixed point. Fig. 5*A*and*B*shows that there is a small excursion beyond the fixed point as the trajectory spirals inward to the new steady-state DIC concentration.Fig. 5.

Next we repeat the same procedure, but with a slightly larger ${\text{CO}}_{2}$ injection rate, $\nu =0.40$. The initial condition remains the same and the new, stable steady-state DIC concentration is only marginally higher than before. However, Fig. 5

*C*and*D*shows that the excursion beyond the fixed point is now roughly an order of magnitude greater than before. Indeed, its size is comparable to the large limit cycle of Fig. 4. The trajectory, however, is not periodic, and the system soon relaxes back to its stable steady state.*SI Appendix*, Fig. S6 shows how the size of such excursions varies with $\nu $. The size is defined by the difference between the maximum values of $w$ and ${w}_{\nu}^{*}$, where ${w}_{\nu}^{*}$ is the DIC fixed point computed from Eq.

**15**for a given value of $\nu $.

*SI Appendix*, Fig. S6 shows that the transition from small to large excursions is sharp yet continuous. Here it occurs near $\nu =0.391$. Beyond this threshold, the size of the excursion no longer grows appreciably with $\nu $.

Dynamical systems that exhibit phenomena similar to those of Fig. 5 and

*SI Appendix*, Fig. S6 are called excitable (21). When a stimulus such as $\nu $ exceeds a threshold, excitable systems produce a transient excitation significantly larger than the stimulus, as seen here. Excitable systems were first introduced to understand the behavior of action potentials (nerve impulses or spiking) in neurons (68, 69). They have since been studied in the context of numerous other problems, including glacial cycles (32) and the response of peatlands to warming (70).Excitations occur near Hopf bifurcations (21). Even though the fixed point is stable, a remnant of the stable limit cycle exists in the phase plane. Here, as illustrated in Fig. 5, when $\nu $ is above a critical threshold ${\nu}_{c}$, the trajectory effectively hops on to a vestige of the limit cycle before settling in to the fixed point. The size and timescale of excitations are therefore comparable to limit cycles. Specifically, their amplitudes are insensitive to the extent to which the threshold parameter exceeds the threshold (

*SI Appendix*, Fig. S6), and their duration is comparable to the period of limit cycles. Unlike limit cycles, however, they do not require that a parameter of the system, such as $\mu $, be increased beyond a point of bifurcation. For the carbon-cycle model studied here, this means that excitations can occur in a stable carbon cycle that is merely subjected to an above-threshold injection of ${\text{CO}}_{2}$ while its constituent parameters—e.g., riverine fluxes into the system, burial fluxes out of it, and the timescale of homeostatic equilibration—remain the same.Fig. 6 shows where excitations occur and at what value of ${\nu}_{c}$, in the plane of ${c}_{x}$ and $\theta $, along with the stability boundary at which the Hopf bifurcation occurs. Fig. 6 also shows the bistable region adjacent to the stability boundary where the subcritical Hopf bifurcation occurs. One sees that excitations occur in a wide region of parameter space adjacent to the bistable region.

Fig. 6.

The relation between excitations and bistability is further clarified in Fig. 7, which provides a 1D view along the ${c}_{x}$ axis. The bistable region occurs at values of ${c}_{x}$ where the stable fixed point and stable limit cycle coexist with an unstable limit cycle. In this region, values of $\nu $ that depress the minimum $w$ of the stable limit cycle below the minimum $w$ of the unstable limit cycle create an immediate jump to the stable limit cycle. At the left boundary of the bistable region, the unstable limit cycle and the stable limit cycle collide via a saddle-node bifurcation of cycles (29). For values of ${c}_{x}$ below this point, only the fixed point is stable. Whereas an above-threshold $\nu $ creates a jump to the stable limit cycle when ${c}_{x}$ is in the bistable region, a similar above-threshold $\nu $ creates an excitation at smaller values of ${c}_{x}$. Fig. 8 illustrates how the unstable limit cycle determines the evolution of phase-plane trajectories in the bistable regime.

*SI Appendix*, Fig. S7 shows that excitations near the saddle-node bifurcation of cycles can create a finite train of pulses before returning to the stable fixed point.Fig. 7.

Fig. 8.

## Discussion

If the carbon cycle is indeed excitable, the great events of the geologic past should reveal signs of such dynamics. Among these events, those associated with mass extinction are of particular interest. It is also of interest to ask how the hypothesis of excitability informs our understanding of the risks of human-induced carbon-cycle disruptions in the future. The following discussion provides some first steps in these directions.

### Characteristic Events.

Because excitations inherit the properties of limit cycles, Eq. where $\alpha $ is a constant. Because $\alpha {j}_{\mathrm{\text{in}}}/{w}^{*}$ is a concentration flux divided by a steady-state concentration, I call it the specific rate of the excitation.

**19**relates not only the amplitude and period of limit cycles via the flux ${j}_{\mathrm{\text{in}}}$, but also the amplitude and duration of excitations. As expected, numerical simulations of the carbon-cycle model show that the excitation size $\mathrm{\Delta}w\simeq {j}_{\mathrm{\text{in}}}\tau $ for large $\tau $, where $\tau $ is the time from the onset of an excitation to its peak (*SI Appendix*, Fig. S8). The geochemical record provides data from which the normalized size $\mathrm{\Delta}w/{w}^{*}$ may be estimated for real events; it also provides estimates of $\tau $ (7). If real events are excitations, then the above reasoning suggests that their size and timescale should approximately satisfy$$\frac{\mathrm{\Delta}w}{{w}^{*}}=\frac{\alpha {j}_{\mathrm{\text{in}}}}{{w}^{*}}\tau ,$$

[20]

Fig. 9 plots $\mathrm{\Delta}w/{w}^{*}$ as a function of $\tau $ for 31 global disruptions of the carbon cycle during the last 542 My (7). Five of the events are associated with major mass extinctions. Variations in $\tau $ are likely related to changes in the timescale of homeostatic damping over geologic time (7). The straight line and its uncertainty (the yellow region) predict the size of the disruption that would occur if it were quantitatively related to the cessation of organic carbon burial and the resulting accumulation of ${\text{CO}}_{2}$ (7). Roughly half of the events occur in this region; they each display surges of carbon that grow with approximately the same specific rate. In addition to being well populated, the yellow region also appears to separate benign events from 4 of the 5 mass extinction events. The disruptions within the yellow region therefore appear to share important characteristics in addition to their common specific rate. Henceforth these disruptions are called characteristic events.

Fig. 9.

I hypothesize that characteristic events represent excitations of the marine carbon cycle. If so, then the size of characteristic events is predicted by Eq.

**20**, just as it is for model excitations. The yellow region in Fig. 9 satisfies Eq.**20**for $\alpha \sim 0.1$ (7). Model excitations scale similarly, but with $\alpha \sim 1$ (*SI Appendix*, Fig. S8). This difference likely reflects mechanisms, such as limitations imposed by organic carbon burial (7), that are not captured by the model. The important point, however, is that for both observed disruptions and model excitations, the amplitude $\mathrm{\Delta}w\propto {j}_{\mathrm{\text{in}}}\tau $.The shared scaling of characteristic events and model excitations with the riverine flux ${j}_{\mathrm{\text{in}}}$ points to a common foundation. The carbon-cycle model is an open, nonlinear chemical reactor driven out of equilibrium by a flux of reactants into and out of the system. So, too, is the real ocean; elemental concentrations are determined by the balance of inputs and outputs while buffering is controlled by heterogeneous carbonate equilibria (71). More generally, instability, bistability, and oscillations can occur in nonlinear chemical systems forced out of equilibrium (72). In such problems, the size of limit cycles—and therefore the size of any excitations—is a property of the system rather than its perturbation (29). The model expressed by Eqs.

**11**and**12**shows how this works in an idealized carbon cycle. The proportionality of event sizes to the driving flux is a robust property of the model that should be equally applicable in more complex settings. The common scaling with ${j}_{\mathrm{\text{in}}}$ also provides strong evidence that characteristic events display an intrinsic nonlinear response to perturbation rather than a proportionate response to an external stress. Because different above-threshold stressors would yield a similar response in an excitable system, the excitation hypothesis also explains why the yellow region in Fig. 9 is relatively well populated. In addition, because the carbon isotopic signals typically return to their initial steady state after each event (e.g., Fig. 1), the steady state apparently remains stable. The combination of these observations is consistent with excitability.Events lying below the yellow region in Fig. 9 occur at relatively slow specific rates and do not appear to share anything in common other than their slowness. They likely represent a quasistatic (10) response to subthreshold changes in the parameters of the carbon cycle, such as changes in burial rates. These minor events therefore probably embody the classic proportionate response. The anomalously fast mass extinction events in the upper region should, however, be related to excitability. To understand how, I first discuss the excitation threshold.

### The Excitation Threshold.

Practically speaking, the amplitude of an excitation should significantly exceed the amplitude of the perturbation required to exceed the threshold (21). In terms of the carbon-cycle model, such a definition immediately provides a conservative upper bound for ${\nu}_{c}$ by requiring that the increase in ${w}^{*}$ due to $\nu $ be smaller than the excitation amplitude $\mathrm{\Delta}w$. Obtaining the former from Eq. Fig. 6 shows that the excitation threshold in the carbon-cycle model is consistent with Eq.

**15**and the latter from Eq.**20**, one then has ${\nu}_{c}\mu <\alpha {j}_{\mathrm{\text{in}}}\tau $ with $\alpha \sim 1$. Similar reasoning applies to the real carbon cycle. Given a damping timescale ${\tau}_{w}$, the DIC fixed point should again increase by roughly ${\nu}_{c}\mu $ at the threshold. The assumption that characteristic events are excitations then provides the same bound on this increase, but with $\alpha \sim 0.1$. In either case, the excitation timescale $\tau $ should be similar to the damping timescale. Recalling $\mu ={j}_{\mathrm{\text{in}}}{\tau}_{w}$, one then finds ${\nu}_{c}<\alpha $ or**21**near the subcritical Hopf bifurcation.The combined output of ${\text{CO}}_{2}$ from modern subaerial and submarine volcanism produces about 0.06 Pg$\cdot \mathrm{C}\cdot {\mathrm{y}}^{-1}$ (73), which is about 15% of the modern riverine flux of $0.41\text{Pg C}\cdot {\text{y}}^{-1}$ (74). If the volcanic flux were doubled, the perturbation expressed by $\nu $ would therefore be about 0.15, equivalent to the predicted upper bound for ${\nu}_{c}$ (Eq.

**22**) and therefore enough, in principle, to initiate an excitation.### Mass Extinctions.

Fig. 9 shows that 4 mass extinction events exhibit a specific rate significantly greater than that of characteristic events, by a factor of 2 to 5 (7). These events can be characterized as excitations driven by an injection rate $\nu $ significantly greater than ${\nu}_{c}$.

The idea follows reasoning similar to that used to bound ${\nu}_{c}$. Excitations at threshold exhibit the characteristic size $\mathrm{\Delta}w\simeq \alpha {j}_{\mathrm{\text{in}}}\tau $ (Eq. In other words, above-threshold excitations are larger than threshold excitations by a factor of about $1+\epsilon {\tau}_{w}/\alpha \tau $. Since the excitation timescale $\tau $ should not appreciably change, the observed rate would similarly increase. Because $\tau \sim {\tau}_{w}$, the main interest lies in $\epsilon /\alpha $, which is on the order of $10\epsilon $ for the real carbon cycle. So, for ${\nu}_{c}\sim 0.1$ (the upper bound of Eq.

**20**). Suppose that ${\text{CO}}_{2}$ injection occurs at rate $\left({\nu}_{c}+\epsilon \right){j}_{\mathrm{\text{in}}}$, where $\epsilon >0$. Eq.**15**then shows that $\mathrm{\Delta}w$ will increase by $\epsilon \mu $ relative to its size at the threshold, resulting in$$\mathrm{\Delta}w\simeq {j}_{\mathrm{\text{in}}}\left(\alpha \tau +\epsilon {\tau}_{w}\right).$$

[23]

**22**), excitation at twice the threshold would increase $\mathrm{\Delta}w$ and its growth rate by roughly a factor of 2.The association of massive flood-basalt volcanism with the end-Permian (23), end-Triassic (22), and end-Cretaceous (24, 25) mass extinctions suggests how this would work. Flood basalts release ${\text{CO}}_{2}$ at a rate of about $3.5\times 1{0}^{-3}\text{Pg C}$/ ${\text{km}}^{3}$ of magma (75). Schoene

*et al.*(24) estimate that up to $11{\text{km}}^{3}\cdot {\text{y}}^{-1}$ of magma erupted from Deccan volcanism over a 10-ky period tens of thousands of years before the end-Cretaceous extinction. The resulting ${\text{CO}}_{2}$ flux is equivalent to about 10% of the modern riverine flux, yielding $\nu \simeq 0.1$, commensurate with the upper bound for ${\nu}_{c}$. Other episodes of massive volcanism (22–25) may exhibit similar pulses (76). Such events may indeed be well above threshold: If ${\nu}_{c}$ were much smaller than its upper bound, then $\epsilon =\nu -{\nu}_{c}\simeq 0.1$, yielding $\alpha /\epsilon \simeq 1$. The size and rate of the event would then be roughly double those of characteristic disruptions.### Excitation of the Modern Carbon Cycle.

As a result of anthropogenic ${\text{CO}}_{2}$ emissions, Earth’s oceans currently absorb about $2.3\hspace{0.17em}\text{Pg C}\cdot {\text{y}}^{-1}$ (77), which is about 5.6 times the riverine carbon influx. The resulting excitation parameter, $\nu =5.6$, is more than 50 times greater than the predicted upper bound for ${\nu}_{c}$. That prediction, however, assumes a continuous, constant injection of ${\text{CO}}_{2}$. If ${\text{CO}}_{2}$ injection occurs instead for a time much less than the damping timescale, its ability to initiate an excitation may be attenuated by the ocean’s homeostatic response.

Whether or not excitation ensues will nevertheless still depend on the injection rate. Suppose that ${\text{CO}}_{2}$ injection occurs over a time ${t}_{i}<{\tau}_{w}$. If excitations follow when $\nu >{\nu}_{c}$ for ${t}_{i}\ge {\tau}_{w}$, then excitations should still occur when ${t}_{i}<{\tau}_{w}$ if the DIC addition $\nu {j}_{\mathrm{\text{in}}}{t}_{i}$ is approximately equal to or greater than the total input resulting from injection at the threshold ${\nu}_{c}$ over 1 damping time. In other words, there is a ${t}_{i}$-dependent threshold ${\nu}_{c}^{\prime}$ that approximately satisfies ${\nu}_{c}^{\prime}{t}_{i}={\nu}_{c}{\tau}_{w}$ for small ${t}_{i}$. Therefore,The scaling like ${t}_{i}^{-1}$ should increase in accuracy as ${t}_{i}\to 0$.

$${\nu}_{c}^{\prime}\left({t}_{i}\right)\simeq {\nu}_{c}\frac{{\tau}_{w}}{{t}_{i}},\text{\hspace{1em}\hspace{1em}}{t}_{i}<{\tau}_{w}.$$

[24]

*SI Appendix*, Fig. S9 shows that numerical solutions of the carbon-cycle model compare well to this prediction. The dependence of the short-time threshold on ${t}_{i}$ is qualitatively similar to the notion of a critical ramping rate (70).Given the real-ocean upper bound ${\nu}_{c}\lesssim 0.1$ (Eq.

**22**), substitution of ${\tau}_{w}=1{0}^{4}$ y and ${t}_{i}=1{0}^{2}$ y into Eq.**23**yields the century-scale upper bound ${\nu}_{c}^{\prime}\lesssim 10$. Mean 21st-century oceanic carbon uptake will likely range from 1.8 to 4.2 $\text{Pg C}\cdot {\text{y}}^{-1}$ (77), which translates to $4.4\le \nu \le 10.2$, similar to the upper bound for ${\nu}_{c}^{\prime}$. With respect to the excitation threshold, the modern perturbation is therefore approximately equivalent to the aforementioned end-Cretaceous pulse; its duration is 2 orders of magnitude shorter but its flux is 2 orders of magnitude greater. Fig. 10 illustrates this correspondence graphically.Fig. 10.

The short-time critical rate ${\nu}_{c}^{\prime}$ may be alternatively expressed as a critical mass ${m}_{c}$. Multiplication of both sides of Eq. where $V=1.35\times 1{0}^{21}\text{kg}$ is the mass of the oceans (2) and the factor of 12 converts mol to g. Ref. 7 derives an equivalent expression from different assumptions. In conjunction with Eq.

**24**by ${j}_{\mathrm{\text{in}}}{t}_{i}$ shows that injection with $\nu >{\nu}_{c}^{\prime}$ over a timescale ${t}_{i}<{\tau}_{w}$ is equivalent to adding a mass of carbon greater than$${m}_{c}\simeq 12{\nu}_{c}{j}_{\mathrm{\text{in}}}{\tau}_{w}V,$$

[25]

**22**, Eq.**25**predicts the upper bound ${m}_{c}\lesssim 405$ Pg C, which falls within the projected uptake of 290 to 540 Pg C from the year 1850 to 2100 (77). Thus, in terms of either 21st-century uptake rates or total postindustrial uptake, the modern marine carbon cycle may be at the precipice of excitation.Finally, there is an important caveat. Since the carbon-cycle model applies only to timescales greater than that of ocean mixing [about 1,000 y (2, 3)], it does not necessarily inform understanding of phenomena at timescales shorter than 1,000 y. However, if the addition of ${\text{CO}}_{2}$ to the oceans were not appreciably damped at submillenial timescales by processes not discussed here, the reasoning outlined above would remain valid at these short timescales.

## Conclusion

Major disruptions of Earth’s carbon cycle are typically interpreted as a proportionate response to an environmental perturbation (8, 16, 20). This paper suggests instead that many of these events represent the nonlinear amplification of processes that operate within the carbon cycle. To illustrate how this could work, this paper constructs and studies a simple dynamical-system model of the marine carbon cycle. The results suggest that disruptions can follow when a stable steady state is perturbed beyond a threshold. A likely source of perturbations is an increased flux of ${\text{CO}}_{2}$ into the oceans. The autocatalytic amplification of model disruptions exhibits a characteristic increase and rate of growth of the ocean’s store of carbon before returning to the original steady state. Such characteristic sizes and rates, independent of the history of external perturbations, are classical features of nonlinear systems (29). The historical record of the real carbon cycle also exhibits similar characteristics (7). These shared properties suggest that the model displays important features of the real carbon cycle.

These findings suggest that a characteristic excitation of upper-ocean ${\text{CO}}_{2}$ levels follows the injection of ${\text{CO}}_{2}$ into the oceans at an above-threshold rate. If injection continues over a time greater than the timescale ${\tau}_{w}$ over which the oceans homeostatically adjust to changes in pH, the threshold’s upper bound is roughly equivalent to the higher estimates of the rate at which ${\text{CO}}_{2}$ degasses from major flood-basalt eruptions associated with mass extinction (24, 76). If injection is instead limited to a shorter duration ${t}_{i}$, the threshold increases by a factor of ${\tau}_{w}/{t}_{i}$. Thus, the relatively slow rate of ${\text{CO}}_{2}$ injection commensurate with massive volcanism at geologic timescales turns out to be roughly equivalent, in terms of its potential to reach the threshold, to the much stronger but briefer perturbation of the modern carbon cycle. Mass extinction events appear to be associated with excitations well above threshold.

These conclusions rely in part on the assumption that the dynamics of a 2D dynamical system represent a subset of the complex behavior exhibited by the real carbon cycle. This assumption does not require that the mechanisms encoded by the terms within the dynamical system be strictly correct. However, it does demand that the model’s qualitative dynamical properties, such as the possible existence of a limit cycle, be realistic. This may not be true. One could argue, for example, that phenomena at submillenial timescales, such as ecological reorganization, will act as strong negative feedbacks in local environments, impeding global autocatalytic self-organization. The present state of the geochemical record makes such notions difficult to rule out observationally. On the other hand, numerical simulation of more detailed carbon-cycle models (44) should allow one to test whether excitations are expressed in such frameworks. Likewise, careful analysis of individual disruptions or specific periods in the geochemical record should allow one to test whether quantitative signatures of excitable systems exist. The possibility of resonance with the fluctuations of external perturbations (21), such as those induced by orbital variations (12), offers one such route. The work reported here identifies why the carbon cycle may be excitable, how excitability may have been expressed in the past, and why an excitation may occur in the future. In a curious twist of geological and biological evolution, dynamical mechanisms underlying environmental upheaval and mass extinction may be similar to those that make neurons spike.

## Acknowledgments

I thank C. Arnscheidt, O. Devauchelle, R. Ferrari, and J. Weitz for helpful discussions and B. Schoene and B. Keller for providing data. This work was supported by NASA Astrobiology Grant NNA13AA90A and NSF Grant EAR-1338810.

## Supporting Information

Appendix (PDF)

- Download
- 1.10 MB

## References

1

R. A. Berner,

*The Phanerozoic Carbon Cycle: CO*(Oxford University Press, New York, 2004)._{2}and O_{2}.2

J. L. Sarmiento, N. Gruber,

*Ocean Biogeochemical Dynamics*(Princeton University Press, Princeton, NJ, 2006).3

S. R. Emerson, J. I. Hedges,

*Chemical Oceanography and the Marine Carbon Cycle*(Cambridge University Press, New York, NY, 2008).4

D. H. Rothman, Earth’s carbon cycle: A mathematical perspective.

*Bull. Am. Math. Soc.***52**, 47–64 (2015).5

O. H. Walliser, Ed.,

*Global Events and Event Stratigraphy in the Phanerozoic*(Springer, Berlin, Germany, 1996).6

S. M. Stanley, Relation of Phanerozoic stable isotope excursions to climate, bacterial metabolism, and major extinctions.

*Proc. Natl. Acad. Sci. U.S.A.***107**, 19185–19189 (2010).7

D. H. Rothman, Thresholds of catastrophe in the Earth system.

*Sci. Adv.***3**, e1700906 (2017).8

L. R. Kump, M. A. Arthur, Interpreting carbon-isotope excursions: Carbonates and organic matter.

*Chem. Geol.***161**, 181–198 (1999).9

J. M. Hayes, H. Strauss, A. J. Kaufman, The abundance of

^{13}C in marine organic matter and isotopic fractionation in the global biogeochemical cycle of carbon during the past 800 Ma.*Chem. Geol.***161**, 103–125 (1999).10

D. H. Rothman, J. M. Hayes, R. E. Summons, Dynamics of the Neoproterozoic carbon cycle.

*Proc. Natl. Acad. Sci. U.S.A.***100**, 8124–8129 (2003).11

P. F. Sexton et al., Eocene global warming events driven by ventilation of oceanic dissolved organic carbon.

*Nature***471**, 349–352 (2011).12

S. K. Turner, P. F. Sexton, C. D. Charles, R. D. Norris, Persistence of carbon release events through the peak of early Eocene global warmth.

*Nat. Geosci.***7**, 748–751 (2014).13

G. R. Dickens, J. R. O’Neil, D. K. Rea, R. M. Owen, Dissociation of oceanic methane hydrate as a cause of the carbon isotope excursion at the end of the Paleocene.

*Paleoceanography***10**, 965–971 (1995).14

L. Alegret, E. Thomas, K. C. Lohmann, End-Cretaceous marine mass extinction not caused by productivity collapse.

*Proc. Natl. Acad. Sci. U.S.A.***109**, 728–732 (2012).15

V. A. Melezhik et al., “The palaeoproterozoic perturbation of the global carbon cycle: The Lomagundi-Jatuli isotopic event” in

*Reading the Archive of Earth’s Oxygenation*, Melezhik V, et al., Eds. (Frontiers in Earth Sciences, Springer, 2013), pp. 1111–1150.16

D. H. Rothman et al., Methanogenic burst in the end-Permian carbon cycle.

*Proc. Natl. Acad. Sci. U.S.A.***111**, 5462–5467 (2014).17

L. R. Kump et al., A weathering hypothesis for glaciation at high atmospheric p${\text{CO}}_{2}$ during the Late Ordovician.

*Palaeogeogr. Palaeoclimatol. Palaeoecol.***152**, 173–187 (1999).18

M. A. Arthur, W. E. Dean, L. M. Pratt, Geochemical and climatic effects of increased marine organic-carbon burial at the Cenomanian/Turonian boundary.

*Nature***335**, 714–717 (1988).19

G. Suan et al., Duration of the Early Toarcian carbon isotope excursion deduced from spectral analysis: Consequence for its possible causes.

*Earth Planet Sci. Lett.***267**, 666–679 (2008).20

Y. Cui, L. R. Kump, A. Ridgwell, Initial assessment of the carbon emission rate and climatic consequences during the end-Permian mass extinction.

*Palaeogeogr. Palaeoclimatol. Palaeoecol.***389**, 128–136 (2013).21

E. M. Izhikevich,

*Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting*(MIT Press, Cambridge, MA, 2007).22

T. J. Blackburn et al., Zircon U-Pb geochronology links the end-Triassic extinction with the central Atlantic magmatic province.

*Science***340**, 941–945 (2013).23

S. D. Burgess, S. A. Bowring, High-precision geochronology confirms voluminous magmatism before, during, and after Earth’s most severe extinction.

*Sci. Adv.***1**, e1500470 (2015).24

B. Schoene et al., U-Pb constraints on pulsed eruption of the Deccan Traps across the end-Cretaceous mass extinction.

*Science***363**, 862–866 (2019).25

C. J. Sprain et al., The eruptive tempo of Deccan volcanism in relation to the Cretaceous-Paleogene boundary.

*Science***363**, 866–870 (2019).26

T. Westerhold et al., Astronomical calibration of the Ypresian timescale: Implications for seafloor spreading rates and the chaotic behavior of the solar system?

*Clim. Past***13**, 1129–1152 (2017).27

U. Riebesell, A. Körtzinger, A. Oschlies, Sensitivities of marine carbon fluxes to ocean change.

*Proc. Natl. Acad. Sci. U.S.A.***106**, 20602–20609 (2009).28

J. Guckenheimer, P. Holmes,

*Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields*(Springer-Verlag, New York, NY, 1983).29

S. Strogatz,

*Nonlinear Dynamics and Chaos*(Addison-Wesley, New York, NY, 1994).30

M. Ghil, S. Childress,

*Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics, Dynamo Theory, and Climate Dynamics*(Springer-Verlag, 1987),**vol. 60**.31

B. Saltzman,

*Dynamical Paleoclimatology: Generalized Theory of Global Climate Change*(Academic Press, 2002).32

M. Crucifix, Oscillators and relaxation phenomena in Pleistocene climate theory.

*Philos. Trans. R. Soc. A***370**, 1140–1165 (2012).33

R. Zeebe, D. Wolf-Gladrow,

*CO*(Elsevier Science, 2001)._{2}in Seawater: Equilibrium, Kinetics, Isotopes34

D. Archer,

*The Long Thaw*(Princeton University Press, Princeton, NJ, 2009).35

D. Archer,

*The Global Carbon Cycle*(Princeton University Press, 2010).36

D. Archer, H. Kheshgi, E. Maier-Reimer, Dynamics of fossil fuel CO

_{2}neutralization by marine CaCO_{3}.*Glob. Biogeochem. Cycles***12**, 259–276 (1998).37

D. Archer, Fate of fossil fuel ${\text{CO}}_{2}$ in geologic time.

*J. Geophys. Res.***110**, C09S05 (2005).38

A. Ridgwell, J. Hargreaves, Regulation of atmospheric ${\text{CO}}_{2}$ by deep-sea sediments in an Earth system model.

*Glob. Biogeochem. Cycles***21**, GB2008 (2007).39

R. E. Zeebe, P. Westbroek, A simple model for the CaCO3 saturation state of the ocean: The “Strangelove,” the “Neritan,” and the “Cretan” Ocean.

*Geochem. Geophys. Geosystems***4**, 1104 (2003).40

S. A. Henson, R. Sanders, E. Madsen, Global patterns in efficiency of particulate organic carbon export and transfer to the deep ocean.

*Glob. Biogeochem. Cycles***26**, GB1028 (2012).41

C. M. Marsay et al., Attenuation of sinking particulate organic carbon flux through the mesopelagic ocean.

*Proc. Natl. Acad. Sci. U.S.A.***112**, 1089–1094 (2015).42

B. B. Cael, M. J. Follows, On the temperature dependence of oceanic export efficiency.

*Geophys. Res. Lett.***43**, 5170–5175 (2016).43

J. Wohlers et al., Changes in biogenic carbon flow in response to sea surface warming.

*Proc. Natl. Acad. Sci U.S.A.***106**, 7067–7072 (2009).44

D. Hülse, S. Arndt, J. D. Wilson, G. Munhoven, A. Ridgwell, Understanding the causes and consequences of past marine carbon cycling variability through models.

*Earth-Sci. Rev.***171**, 349–382 (2017).45

U. Riebesell et al., Reduced calcification of marine plankton in response to increased atmospheric CO

_{2}.*Nature***407**, 364–367 (2000).46

A. Sciandra et al., Response of coccolithophorid Emiliania huxleyi to elevated partial pressure of ${\text{CO}}_{2}$ under nitrogen limitation.

*Mar. Ecol. Prog. Ser.***261**, 111–122 (2003).47

B. Delille et al., Response of primary production and calcification to changes of p${\text{CO}}_{2}$ during experimental blooms of the coccolithophorid Emiliania huxleyi.

*Glob. Biogeochem. Cycles***19**, GB2023 (2005).48

Y. Feng et al., Interactive effects of increased p${\text{CO}}_{2}$, temperature and irradiance on the marine coccolithophore Emiliania huxleyi (Prymnesiophyceae).

*Eur. J. Phycol.***43**, 87–98 (2008).49

R. Francois, S. Honjo, R. Krishfield, S. Manganini, Factors controlling the flux of organic carbon to the bathypelagic zone of the ocean.

*Glob. Biogeochem. Cycles***16**, 34-1–34-20 (2002).50

R. A. Armstrong, C. Lee, J. I. Hedges, S. Honjo, S. G. Wakeham, A new, mechanistic model for organic carbon fluxes in the ocean based on the quantitative association of POC with ballast materials.

*Deep-Sea Res.***49**, 219–236 (2002).51

C. Klaas, D. E. Archer, Association of sinking organic matter with various types of mineral ballast in the deep sea: Implications for the rain ratio.

*Glob. Biogeochem. Cycles***16**, 63-1–63-14 (2002).52

S. Barker, J. A. Higgins, H. Elderfield, The future of the carbon cycle: Review, calcification response, ballast and feedback on atmospheric CO

_{2}.*Philos. Trans. R. Soc. Lond. A***361**, 1977–1999 (2003).53

M. Hofmann, H. J. Schellnhuber, Oceanic acidification affects marine carbon pump and triggers extended marine oxygen holes.

*Proc. Natl. Acad. Sci. U.S.A.***106**, 3017–3022 (2009).54

U. Passow, Switching perspectives: Do mineral fluxes determine particulate organic carbon fluxes or vice versa?

*Geochem. Geophys. Geosystems***5**, Q04002 (2004).55

P. J. Lam, S. C. Doney, J. K. B. Bishop, The dynamic ocean biological pump: Insights from a global compilation of particulate organic carbon, ${\text{CaCO}}_{3}$, and opal concentration profiles from the mesopelagic.

*Glob. Biogeochem. Cycles***25**, GB3009 (2011).56

J. Wilson, S. Barker, A. Ridgwell, Assessment of the spatial variability in particulate organic matter and mineral sinking fluxes in the ocean interior: Implications for the ballast hypothesis.

*Glob. Biogeochem. Cycles***26**, GB4011 (2012).57

F. A. C. Le Moigne et al., On the proportion of ballast versus non-ballast associated carbon export in the surface ocean.

*Geophys. Res. Lett.***39**, L15610 (2012).58

F. A. C. Le Moigne, K. Pabortsava, C. L. J. Marcinko, P. Martin, R. J. Sanders, Where is mineral ballast important for surface export of particulate organic carbon in the ocean?

*Geophys. Res. Lett.***41**, 8460–8468 (2014).59

T. DeVries, T. Weber, The export and fate of organic matter in the ocean: New constraints from combining satellite and oceanographic tracer observations.

*Glob. Biogeochem. Cycles***31**, 535–555 (2017).60

C. Laufkötter, J. G. John, C. A. Stock, J. P. Dunne, Temperature and oxygen dependence of the remineralization of organic matter.

*Glob. Biogeochem. Cycles***31**, 1038–1050 (2017).61

J. A. Cram et al., The role of particle size, ballast, temperature, and oxygen in the sinking flux to the deep sea.

*Glob. Biogeochem. Cycles***32**, 858–876 (2018).62

J. P. Gattuso, L. Hansson,

*Ocean Acidification*(Oxford University Press, 2011).63

P. F. Hoffman, D. P. Schrag, The snowball Earth hypothesis: Testing the limits of global change.

*Terra Nova***14**, 129–155 (2002).64

I. C. Handoh, T. M. Lenton, Periodic mid-Cretaceous oceanic anoxic events linked by oscillations of the phosphorus and oxygen biogeochemical cycles.

*Glob. Biogeochem. Cycles***17**, 1092 (2003).65

A. C. Maloof, D. P. Schrag, J. L. Crowley, S. A. Bowring, An expanded record of Early Cambrian carbon cycling from the Anti-Atlas Margin, Morocco.

*Can. J. Earth Sci.***42**, 2195–2216 (2005).66

J. L. Payne et al., Large perturbations of the carbon cycle during recovery from the end-Permian extinction.

*Science***305**, 506–509 (2004).67

A. Bachan et al., A model for the decrease in amplitude of carbon isotope excursions across the phanerozoic.

*Am. J. Sci.***317**, 641–676 (2017).68

A. L. Hodgkin, A. F. Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve.

*J. Physiol.***117**, 500–544 (1952).69

R. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane.

*Biophys. J.***1**, 445–466 (1961).70

S. Wieczorek, P. Ashwin, C. M. Luke, P. M. Cox, Excitability in ramped systems: The compost-bomb instability.

*Proc. R. Soc. A***467**, 1243–1269 (2011).71

R. E. McDuff, F. M. M. Morel, The geochemical control of seawater (Sillen revisited).

*Environ. Sci. Technol.***14**, 1182–1186 (1980).72

J. Ross, Temporal and spatial structures in chemical instabilities.

*Berichte der Bunsengesellschaft für Physikalische Chem.***80**, 1112–1125 (1976).73

P. Huybers, C. Langmuir, Feedback between deglaciation, volcanism, and atmospheric ${\mathrm{c}\mathrm{o}}_{2}$.

*Earth Planet Sci. Lett.***286**, 479–491 (2009).74

W. J. Cai et al., A comparative overview of weathering intensity and ${\text{HCO}}_{3}^{-}$ flux in the world’s major rivers with emphasis on the Changjiang, Huanghe, Zhujiang (Pearl) and Mississippi rivers.

*Cont. Shelf Res.***28**, 1538–1549 (2008).75

S. Self, T. Thordarson, M. Widdowson, Gas fluxes from flood basalt eruptions.

*Elements***1**, 283–287 (2005).76

V. Courtillot, F. Fluteau, “A review of the embedded time scales of flood basalt volcanism with special emphasis on dramatically short magmatic pulses” in Volcanism, Impacts, and Mass Extinctions: Causes and Effects, G. Keller, A. Kerr, Eds. (Geological Society of America Special Paper, 2014), vol. 505.

77

P. Ciais et al., “Carbon and other biogeochemical cycles” in

*Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change*, T. F. Stocker et al., Eds. (Cambridge University Press, Cambridge, UK, 2013), pp. 465–570.## Information & Authors

### Information

#### Published in

#### Classifications

#### Copyright

© 2019. Published under the PNAS license.

#### Submission history

**Published online**: July 8, 2019

**Published in issue**: July 23, 2019

#### Keywords

#### Acknowledgments

I thank C. Arnscheidt, O. Devauchelle, R. Ferrari, and J. Weitz for helpful discussions and B. Schoene and B. Keller for providing data. This work was supported by NASA Astrobiology Grant NNA13AA90A and NSF Grant EAR-1338810.

#### Notes

This article is a PNAS Direct Submission.

### Authors

#### Competing Interests

The author declares no conflict of interest.

## Metrics & Citations

### Metrics

#### Citation statements

#### Altmetrics

### Citations

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

#### Cited by

Loading...

## View Options

### View options

#### PDF format

Download this article as a PDF file

DOWNLOAD PDF### Get Access

#### Login options

Check if you have access through your login credentials or your institution to get full access on this article.

Personal login Institutional Login#### Recommend to a librarian

Recommend PNAS to a Librarian#### Purchase options

Purchase this article to get full access to it.