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)
July 8, 2019
116 (30) 14813-14822


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.


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 CO2 enters the oceans at a flux that exceeds a threshold. The threshold depends on the duration of the injection. For injections lasting a time ti10,000 y in the modern carbon cycle, the threshold flux is constant; for smaller ti, the threshold scales like ti1. Consequently the unusually strong but geologically brief duration of modern anthropogenic oceanic CO2 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.


Earth’s carbon cycle is a loop between photosynthesis, which converts carbon dioxide (CO2) to organic carbon, and respiration, which converts organic carbon back to CO2 (14). The cycle has undergone many disruptions throughout Earth’s history (57). 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 (810). 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.
Fluctuations of the isotopic composition of carbonate carbon (δ13C) during the Eocene period, about 54 Ma (26). Time advances to the right and is given with respect to the minimum of the first abrupt downswing, an event known as Eocene Thermal Maximum 2 (ETM2) or H1. The second event, about 100 ky later, is called H2. The timescale is derived from astrochronology (26).
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 CO2. 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 CO2 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 (2225).
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 (2832).
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 (HCO3) and carbonate ions (CO32); the factor of 2 attached to the CO32 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. 3 specify the equilibrium of the carbonate system. There are 6 unknowns: Not only a and w, but also the concentrations of CO2, CO32, HCO3, and H+. Specification of any 2 of these quantities therefore provides the other 4 via equations (13, 33).

General Formulation.

I consider a well-mixed open system in which dissolved calcium carbonate (CaCO3) 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 reads
where dots represent time derivatives. The concentrations a and w are expressed in units of μmolkg1. The factor of 2 arises because 2 molar equivalents of alkalinity are associated with each mole of CaCO3 (33). Expressing the model in terms of the conserved quantities a and w rather than, e.g., the concentrations of CO2 and CO32 makes it possible to temporarily forestall a discussion of the changing chemical equilibrium.
Fig. 2.
Schematic diagram of the model (not drawn to scale). Left and Right panels represent, respectively, the evolution of total alkalinity, a, and total dissolved inorganic carbon (DIC), w. The wavy line represents the air–sea interface, the upper horizontal thick line divides the shallow ocean from the deep sea, and the lower horizontal dashed line represents the sediment–seawater interface. Concentration fluxes are indicated by unidirectional arrows. The sediment–seawater interface is dashed to indicate that there is no dynamical distinction between the deep sea and sediments.
The flux jin represents the rate at which dissolved CaCO3 is added to the oceans by rivers in terms of an equivalent change in the whole-ocean concentration of dissolved CaCO3, expressed in units of μmolkg1y1. The flux B0 denotes the rate at which CaCO3 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 CaCO3 exits the shallow layer without returning. The dimensionless term ν determines the strength of an external source of CO2, such as volcanic emissions.
The flux R0 represents the extent to which upper-ocean production of CO2 by respiration exceeds a theoretical baseline. The small change in alkalinity associated with respiration (33) is neglected. When ν=0, the baseline case R=0 is associated with a steady state in which w=w0 and the riverine influx of dissolved CaCO3 is precisely balanced by burial (i.e., B=jin). 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 w0. The model assumes that this increase is resisted linearly, so that if the baseline R=0 were restored, w would return to w0 at the timescale τw. At long timescales (102 y) in the real ocean, the processes responsible for such dissipation include the riverine influx of dissolved CaCO3 and the dissolution and precipitation of CaCO3 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 jin 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 ww0. The homeostat’s dominant characteristic timescale, denoted here by τw, is about 10 ky (3638) 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 CaCO3 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, CaCO3 dissolves; when it is supersaturated, CaCO3 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 by
where the carbonate ion concentration c(a,w) derives from the equilibrium of the carbonate system (33) and
is the sigmoidal function shown in SI Appendix, Fig. S1. The exponent γ parameterizes the sharpness of the transition, cp is the transitional carbonate ion concentration where s=1/2, and bjin 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 γ may also appear unwelcome; however, this paper’s conclusions require only that the function s be sigmoidal. SI Appendix, Fig. S1 suggests that γ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, 4042). Possible consequences include an upward shift of respiration and a positive feedback for shallow-ocean and atmospheric CO2 levels (27, 41, 43, 44). Alternatively, carbonate production by pelagic calcifiers may decrease at lower pH (4548). The “ballast hypothesis” (4951) 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 (5458) nor the temperature (42, 5961) 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 CO32 lead to less biogenic carbonate production and therefore less ballast, less associated organic carbon export, and more upper-ocean respiration (27, 4553). 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 (SI Appendix, Fig. S1), but with a possibly different amplitude and transitional concentration. I therefore write
where s¯=1s, θjin is the maximum value of R, and cx is the “crossover” concentration that marks the halfway point in the transition from high to low R. As for the burial function of Eq. 6, the sharpness of the transition parameterized by γ 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 ȧ to ċ.

SI Appendix shows that the rate of change of the CO32 concentration can be approximated by
Here the “buffer function”
approximates the magnitude of the partial derivatives c/a and c/w via specific values of f0, β, and cf related to the equilibrium of the carbonate system (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
where time is nondimensionalized by tt/τw and
is a characteristic concentration. The system of Eqs. 11 and 12 is henceforth called the carbon-cycle model. Its parameters are listed in SI Appendix, Table S1. Of these, only μ, b, θ, cx, cp, and ν 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 ċ=ẇ=0. Under this condition, the addition of Eqs. 11 and 12 yields bs(c,cp)=1. Substituting Eq. 7 and solving for c, one then finds the unique steady-state carbonate ion concentration
which exists only for b>1. Substitution of c* for c in Eq. 12 then provides the steady-state DIC concentration
Note that the untransformed model (Eqs. 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 Tr and determinant Δ of the Jacobian matrix of Eqs. 11 and 12 evaluated at the fixed point. These calculations yield
The fixed point is unstable when Tr>0. Fig. 3 provides stability boundaries in the planes of μ and b, μ and cx, and μ and θ. Instability is favored for large μ=jinτw because greater rates of riverine input and longer damping times provide for faster-growing fluctuations. For some combinations of parameters, however, homeostatic damping always maintains stability.
Fig. 3.
(A–C) Stability in the planes of μ and b (A), μ and cx (B), and μ and θ (C). Values of fixed parameters are given in SI Appendix, Table S1. Crossing the red and blue boundaries results in supercritical and subcritical Hopf bifurcations, respectively. The fixed point is stable only where indicated. However, the region of stability of the limit cycle extends into the region of the stable fixed point adjacent to the subcritical bifurcation. Fig. 6 maps the bistable region in the plane of cx and θ.

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(t) and w(t) on the limit cycle, along with the time evolution of quantities derived from these 2 quantities.
Fig. 4.
A stable limit cycle in the phase plane of the CO32 concentration c and the DIC concentration w. Values of parameters are given in SI Appendix, Table S1. The dotted curves represent the nullclines ċ=0 and ẇ=0.
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 ẇ=0. Here, the ballast feedback begins to increase DIC by the addition of CO2, which acts to consume CO32 ions via the reaction of SI Appendix, Eq. S8. However, the relative rate at which the CO32 concentration decreases compared with the rate at which DIC increases becomes progressively slower as the ocean becomes less buffered at lower CO32 concentrations. Eventually the CO32 concentration reaches a minimum where the loss of CO32 via burial and the CO2 feedback is balanced by the homeostatic response, and it crosses the nullcline ċ=0. However, DIC continues to rise until it is overcome by the homeostatic feedback, reaching a maximum where it crosses the nullcline ẇ=0. The homeostatic feedback then acts to decrease w while increasing c. The CO32 concentration reaches its maximum when carbonate burial becomes stronger than the addition of dissolved CaCO3 from rivers and the weakened homeostatic feedback. A rapid upward jump of the burial rate then acts to reduce both the CO32 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 CO2 (SI Appendix, Fig. S4C) and a reduction of pH (SI Appendix, Fig. S4E). 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 CO2 levels (34). Finally, the surfeit of dissolved CaCO3 must eventually be buried. In the limit cycle pictured here, carbonate burial occurs as a sharp pulse (SI Appendix, Fig. S4H) 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 (6467). 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 Δw=wmaxw*, where wmax is the maximum value of w in the cycle; one similarly has Δc=cmaxc*. SI Appendix shows that Δw and Δc generally increase with μ and θ but decrease with b. The limit cycle’s size grows roughly linearly with μ because μ sets the rate of all processes other than homeostatic damping. Because rates and size scale similarly with μ, the dimensionless period T of all but the smallest limit cycles should be effectively independent of μ. These observations imply that
where A=(ΔwΔc)1/2 is a measure of the amplitude of the limit cycle in units of μmolkg1 and T01 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, 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 μ (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 τ=Tτw/2, τ0=T0τw/2, and recalling μ=jinτw, one finds
This relationship, which essentially follows from dimensional analysis, shows that, for ττ0, the limit cycle’s size is roughly given by the riverine influx jin integrated over the time required for the system to be forced maximally out of its steady state.
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.


Now suppose the system is in a stable steady state with no injection of CO2 (ν=0). At time t=0 we set ν=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 CO2 injection increases w* by νμ, but c* is unchanged. The new fixed point retains its stability, because stability is independent of ν (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.
Evolution following perturbations below and above the excitations threshold. (A–D) Phase-plane representation and the time series w(t) for a subthreshold perturbation with ν=0.35 (A and B) and an above-threshold perturbation with ν=0.40 (C and D). The initial condition and all other parameters are the same in each simulation. The crossover CO32 concentration cx=55μmolkg1; other parameters are given in SI Appendix, Table S1. The larger value of ν in the above-threshold case increases the DIC fixed point w* by 0.05μ=12.5μmolkg1 in C compared with A, which is nearly imperceptible in the plots.
Next we repeat the same procedure, but with a slightly larger CO2 injection rate, ν=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 ν. The size is defined by the difference between the maximum values of w and wν*, where wν* is the DIC fixed point computed from Eq. 15 for a given value of ν. SI Appendix, Fig. S6 shows that the transition from small to large excursions is sharp yet continuous. Here it occurs near ν=0.391. Beyond this threshold, the size of the excursion no longer grows appreciably with ν.
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 ν 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 ν is above a critical threshold ν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 μ, 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 CO2 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 νc, in the plane of cx and θ, 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.
Contours of νc (dotted) in the region of excitations (pale yellow) and jumps to the stable limit cycle in the region of bistability (pale blue) in the plane of cx and θ. The blue and red portions of the stability boundary correspond to a subcritical and a supercritical Hopf bifurcation, respectively. Fixed parameters are specified in SI Appendix, Table S1. Excitations or jumps are defined to occur if there exists a ν for which dΔw/dν>103, where Δw=wmaxw*. The threshold νc is the value of ν where that derivative is greatest.
The relation between excitations and bistability is further clarified in Fig. 7, which provides a 1D view along the cx axis. The bistable region occurs at values of cx where the stable fixed point and stable limit cycle coexist with an unstable limit cycle. In this region, values of ν 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 cx below this point, only the fixed point is stable. Whereas an above-threshold ν creates a jump to the stable limit cycle when cx is in the bistable region, a similar above-threshold ν creates an excitation at smaller values of cx. 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.
Bifurcation diagram for the parameter cx, the crossover CO32 concentration for the ballast feedback. The solid and dotted black line respectively represents the stable and unstable fixed point w*. The solid blue line indicates the maximum and minimum values of w in the stable limit cycle and the dashed red line represents the same extremes of the unstable limit cycle. The subcritical Hopf bifurcation occurs where the radius of the unstable limit cycle goes to zero, at cx=62.61μmolkg1. The saddle-node bifurcation of cycles occurs where the unstable and stable limit cycles collide, at cx=55.89μmolkg1. Excitations occur at smaller values of cx, to the left of the arrowhead. Fixed parameters are specified in SI Appendix, Table S1.
Fig. 8.
Phase-space trajectories in the bistable regime. Here cx=58μmolkg1 and all other parameters are the same as in Fig. 5 C and D. The red dashed limit cycle is unstable; thus trajectories initialized inside it return to the (stable) fixed point. Trajectories initialized outside the unstable limit cycle evolve to the stable limit cycle.


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. 19 relates not only the amplitude and period of limit cycles via the flux jin, but also the amplitude and duration of excitations. As expected, numerical simulations of the carbon-cycle model show that the excitation size Δwjinτ for large τ, where τ 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 Δw/w* may be estimated for real events; it also provides estimates of τ (7). If real events are excitations, then the above reasoning suggests that their size and timescale should approximately satisfy
where α is a constant. Because αjin/w* is a concentration flux divided by a steady-state concentration, I call it the specific rate of the excitation.
Fig. 9 plots Δw/w* as a function of τ 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 τ 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 CO2 (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.
The relationship between the relative size and duration of 31 disruptions of the global carbon cycle during the last 542 My (7). The relative size Δw/w* is obtained from changes in the isotopic composition of carbonate carbon that occur in time series similar to that of Fig. 1. The yellow region contains characteristic events; these events satisfy Eq. 20 with α0.1. The events in the pale red region, which include 4 of the 5 mass extinction events, grow faster than characteristic events. The light blue region contains minor events with relatively slow growth rates. The labeled events are associated with the end-Cretaceous (KT), end-Triassic (TJ), end-Permian (PT), end-Ordovician (Ord), and Frasnian–Famennian (FF) mass extinctions. Ref. 7 provides descriptions of each event and a discussion of uncertainties.
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 α0.1 (7). Model excitations scale similarly, but with α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 Δwjinτ.
The shared scaling of characteristic events and model excitations with the riverine flux jin 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 jin 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 νc by requiring that the increase in w* due to ν be smaller than the excitation amplitude Δw. Obtaining the former from Eq. 15 and the latter from Eq. 20, one then has νcμ<αjinτ with α1. Similar reasoning applies to the real carbon cycle. Given a damping timescale τw, the DIC fixed point should again increase by roughly νcμ at the threshold. The assumption that characteristic events are excitations then provides the same bound on this increase, but with α0.1. In either case, the excitation timescale τ should be similar to the damping timescale. Recalling μ=jinτw, one then finds νc<α or
Fig. 6 shows that the excitation threshold in the carbon-cycle model is consistent with Eq. 21 near the subcritical Hopf bifurcation.
The combined output of CO2 from modern subaerial and submarine volcanism produces about 0.06 PgCy1 (73), which is about 15% of the modern riverine flux of 0.41Pg Cy1 (74). If the volcanic flux were doubled, the perturbation expressed by ν would therefore be about 0.15, equivalent to the predicted upper bound for ν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 ν significantly greater than νc.
The idea follows reasoning similar to that used to bound νc. Excitations at threshold exhibit the characteristic size Δwαjinτ (Eq. 20). Suppose that CO2 injection occurs at rate (νc+ε)jin, where ε>0. Eq. 15 then shows that Δw will increase by εμ relative to its size at the threshold, resulting in
In other words, above-threshold excitations are larger than threshold excitations by a factor of about 1+ετw/ατ. Since the excitation timescale τ should not appreciably change, the observed rate would similarly increase. Because ττw, the main interest lies in ε/α, which is on the order of 10ε for the real carbon cycle. So, for νc0.1 (the upper bound of Eq. 22), excitation at twice the threshold would increase Δ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 CO2 at a rate of about 3.5×103Pg C/ km3 of magma (75). Schoene et al. (24) estimate that up to 11km3y1 of magma erupted from Deccan volcanism over a 10-ky period tens of thousands of years before the end-Cretaceous extinction. The resulting CO2 flux is equivalent to about 10% of the modern riverine flux, yielding ν0.1, commensurate with the upper bound for νc. Other episodes of massive volcanism (2225) may exhibit similar pulses (76). Such events may indeed be well above threshold: If νc were much smaller than its upper bound, then ε=ννc0.1, yielding α/ε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 CO2 emissions, Earth’s oceans currently absorb about 2.3Pg Cy1 (77), which is about 5.6 times the riverine carbon influx. The resulting excitation parameter, ν=5.6, is more than 50 times greater than the predicted upper bound for νc. That prediction, however, assumes a continuous, constant injection of CO2. If CO2 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 CO2 injection occurs over a time ti<τw. If excitations follow when ν>νc for tiτw, then excitations should still occur when ti<τw if the DIC addition νjinti is approximately equal to or greater than the total input resulting from injection at the threshold νc over 1 damping time. In other words, there is a ti-dependent threshold νc that approximately satisfies νcti=νcτw for small ti. Therefore,
νc(ti)νcτwti,  ti<τw.
The scaling like ti1 should increase in accuracy as ti0. 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 ti is qualitatively similar to the notion of a critical ramping rate (70).
Given the real-ocean upper bound νc0.1 (Eq. 22), substitution of τw=104 y and ti=102 y into Eq. 23 yields the century-scale upper bound νc10. Mean 21st-century oceanic carbon uptake will likely range from 1.8 to 4.2 Pg Cy1 (77), which translates to 4.4ν10.2, similar to the upper bound for νc. 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.
Equivalence, with respect to the excitation threshold, of the perturbations of the modern and end-Cretaceous carbon cycles, expressed in terms of the dimensionless CO2 injection rate ν. The straight line labeled νc is the threshold’s upper bound predicted by Eq. 24, assuming the modern damping timescale τw=104 y; the segment labeled νc (Eq. 22) provides the upper bound for times ti>104 y. The circles represent projected 21st-century ocean carbon uptake rates for a range of plausible scenarios (77); the vertical line indicates their uncertainty. The symbols labeled KT¯ provide the median and upper limit of the 68% confidence interval estimated for the period of peak Deccan volcanism tens of thousands of years before the end-Cretaceous extinction (24). Because the excitation threshold scales like 1/ti, both the modern and end-Cretaceous perturbations potentially lie near the threshold’s upper bound despite the 2 orders of magnitude that separate their rates.
The short-time critical rate νc may be alternatively expressed as a critical mass mc. Multiplication of both sides of Eq. 24 by jinti shows that injection with ν>νc over a timescale ti<τw is equivalent to adding a mass of carbon greater than
where V=1.35×1021kg 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. 22, Eq. 25 predicts the upper bound mc405 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 CO2 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.


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 CO2 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 CO2 levels follows the injection of CO2 into the oceans at an above-threshold rate. If injection continues over a time greater than the timescale τ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 CO2 degasses from major flood-basalt eruptions associated with mass extinction (24, 76). If injection is instead limited to a shorter duration ti, the threshold increases by a factor of τw/ti. Thus, the relatively slow rate of CO2 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.


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)


R. A. Berner, The Phanerozoic Carbon Cycle: CO2 and O2. (Oxford University Press, New York, 2004).
J. L. Sarmiento, N. Gruber, Ocean Biogeochemical Dynamics (Princeton University Press, Princeton, NJ, 2006).
S. R. Emerson, J. I. Hedges, Chemical Oceanography and the Marine Carbon Cycle (Cambridge University Press, New York, NY, 2008).
D. H. Rothman, Earth’s carbon cycle: A mathematical perspective. Bull. Am. Math. Soc. 52, 47–64 (2015).
O. H. Walliser, Ed., Global Events and Event Stratigraphy in the Phanerozoic (Springer, Berlin, Germany, 1996).
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).
D. H. Rothman, Thresholds of catastrophe in the Earth system. Sci. Adv. 3, e1700906 (2017).
L. R. Kump, M. A. Arthur, Interpreting carbon-isotope excursions: Carbonates and organic matter. Chem. Geol. 161, 181–198 (1999).
J. M. Hayes, H. Strauss, A. J. Kaufman, The abundance of 13C 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).
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).
P. F. Sexton et al., Eocene global warming events driven by ventilation of oceanic dissolved organic carbon. Nature 471, 349–352 (2011).
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).
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).
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).
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.
D. H. Rothman et al., Methanogenic burst in the end-Permian carbon cycle. Proc. Natl. Acad. Sci. U.S.A. 111, 5462–5467 (2014).
L. R. Kump et al., A weathering hypothesis for glaciation at high atmospheric pCO2 during the Late Ordovician. Palaeogeogr. Palaeoclimatol. Palaeoecol. 152, 173–187 (1999).
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).
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).
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).
E. M. Izhikevich, Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting (MIT Press, Cambridge, MA, 2007).
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).
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).
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).
C. J. Sprain et al., The eruptive tempo of Deccan volcanism in relation to the Cretaceous-Paleogene boundary. Science 363, 866–870 (2019).
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).
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).
J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields (Springer-Verlag, New York, NY, 1983).
S. Strogatz, Nonlinear Dynamics and Chaos (Addison-Wesley, New York, NY, 1994).
M. Ghil, S. Childress, Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics, Dynamo Theory, and Climate Dynamics (Springer-Verlag, 1987), vol. 60.
B. Saltzman, Dynamical Paleoclimatology: Generalized Theory of Global Climate Change (Academic Press, 2002).
M. Crucifix, Oscillators and relaxation phenomena in Pleistocene climate theory. Philos. Trans. R. Soc. A 370, 1140–1165 (2012).
R. Zeebe, D. Wolf-Gladrow, CO2 in Seawater: Equilibrium, Kinetics, Isotopes (Elsevier Science, 2001).
D. Archer, The Long Thaw (Princeton University Press, Princeton, NJ, 2009).
D. Archer, The Global Carbon Cycle (Princeton University Press, 2010).
D. Archer, H. Kheshgi, E. Maier-Reimer, Dynamics of fossil fuel CO2 neutralization by marine CaCO3. Glob. Biogeochem. Cycles 12, 259–276 (1998).
D. Archer, Fate of fossil fuel CO2 in geologic time. J. Geophys. Res. 110, C09S05 (2005).
A. Ridgwell, J. Hargreaves, Regulation of atmospheric CO2 by deep-sea sediments in an Earth system model. Glob. Biogeochem. Cycles 21, GB2008 (2007).
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).
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).
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).
B. B. Cael, M. J. Follows, On the temperature dependence of oceanic export efficiency. Geophys. Res. Lett. 43, 5170–5175 (2016).
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).
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).
U. Riebesell et al., Reduced calcification of marine plankton in response to increased atmospheric CO2. Nature 407, 364–367 (2000).
A. Sciandra et al., Response of coccolithophorid Emiliania huxleyi to elevated partial pressure of CO2 under nitrogen limitation. Mar. Ecol. Prog. Ser. 261, 111–122 (2003).
B. Delille et al., Response of primary production and calcification to changes of pCO2 during experimental blooms of the coccolithophorid Emiliania huxleyi. Glob. Biogeochem. Cycles 19, GB2023 (2005).
Y. Feng et al., Interactive effects of increased pCO2, temperature and irradiance on the marine coccolithophore Emiliania huxleyi (Prymnesiophyceae). Eur. J. Phycol. 43, 87–98 (2008).
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).
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).
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).
S. Barker, J. A. Higgins, H. Elderfield, The future of the carbon cycle: Review, calcification response, ballast and feedback on atmospheric CO2. Philos. Trans. R. Soc. Lond. A 361, 1977–1999 (2003).
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).
U. Passow, Switching perspectives: Do mineral fluxes determine particulate organic carbon fluxes or vice versa? Geochem. Geophys. Geosystems 5, Q04002 (2004).
P. J. Lam, S. C. Doney, J. K. B. Bishop, The dynamic ocean biological pump: Insights from a global compilation of particulate organic carbon, CaCO3, and opal concentration profiles from the mesopelagic. Glob. Biogeochem. Cycles 25, GB3009 (2011).
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).
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).
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).
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).
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).
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).
J. P. Gattuso, L. Hansson, Ocean Acidification (Oxford University Press, 2011).
P. F. Hoffman, D. P. Schrag, The snowball Earth hypothesis: Testing the limits of global change. Terra Nova 14, 129–155 (2002).
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).
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).
J. L. Payne et al., Large perturbations of the carbon cycle during recovery from the end-Permian extinction. Science 305, 506–509 (2004).
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).
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).
R. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane. Biophys. J. 1, 445–466 (1961).
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).
R. E. McDuff, F. M. M. Morel, The geochemical control of seawater (Sillen revisited). Environ. Sci. Technol. 14, 1182–1186 (1980).
J. Ross, Temporal and spatial structures in chemical instabilities. Berichte der Bunsengesellschaft für Physikalische Chem. 80, 1112–1125 (1976).
P. Huybers, C. Langmuir, Feedback between deglaciation, volcanism, and atmospheric co2. Earth Planet Sci. Lett. 286, 479–491 (2009).
W. J. Cai et al., A comparative overview of weathering intensity and HCO3 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).
S. Self, T. Thordarson, M. Widdowson, Gas fluxes from flood basalt eruptions. Elements 1, 283–287 (2005).
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.
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


Published in

Go to Proceedings of the National Academy of Sciences
Go to Proceedings of the National Academy of Sciences
Proceedings of the National Academy of Sciences
Vol. 116 | No. 30
July 23, 2019
PubMed: 31285325


Submission history

Published online: July 8, 2019
Published in issue: July 23, 2019


  1. carbon cycle
  2. mass extinctions
  3. excitable systems
  4. dynamical systems
  5. carbon isotopic events


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.


This article is a PNAS Direct Submission.



Daniel H. Rothman1 [email protected]
Lorenz Center, Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139


Author contributions: D.H.R. designed research, performed research, analyzed data, and wrote the paper.

Competing Interests

The author declares no conflict of interest.

Metrics & Citations


Note: The article usage is presented with a three- to four-day delay and will update daily once available. Due to ths delay, usage data will not appear immediately following publication. Citation information is sourced from Crossref Cited-by service.

Citation statements



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


    View Options

    View options

    PDF format

    Download this article as a PDF file


    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.

    Single Article Purchase

    Characteristic disruptions of an excitable carbon cycle
    Proceedings of the National Academy of Sciences
    • Vol. 116
    • No. 30
    • pp. 14783-15309







    Share article link

    Share on social media