SIAM News Blog
SIAM News

# Carbon Cycle Catastrophes: A Dynamical Systems Perspective

Earth’s carbon cycles between photosynthesis—which converts carbon dioxide (CO2) to organic carbon—and respiration, the metabolic processes that transform organic carbon to CO2 [1, 3]. Since pre-industrial times, human activities have caused atmospheric CO2 levels to rise by nearly 50 percent, while roughly half as much CO2 has entered the oceans. While there is considerable interest in predicting climatic responses to these changes, a larger question concerns consequences for the Earth system as a whole — for interactions of life and the environment, in addition to climate. The geologic record indicates that each of the five great mass extinction events of the last 540 million years is associated with significant disruptions of the carbon cycle, as are periods of environmental change unrelated to mass extinction. Although the causes of these events remain ambiguous, the disruptions themselves likely represent a qualitative change in the dynamics of the Earth system. Here I demonstrate how dynamical systems theory helps researchers understand these events.

Figure 1. Disruptions of the carbon cycle. 1a. Time series of the isotopic composition of carbonate carbon ($$\delta ^{13}\mbox{C}$$, the relative enrichment of $$^{13}\mbox{C}$$ compared to $$^{12}\mbox{C}$$, expressed as “per mil”) during the Eocene period, about 54 million years ago. The two abrupt downswings correspond to increases $$\Delta m$$ in the mass $$m$$ of dissolved inorganic carbon in the oceans. 1b. The relative size and duration of 31 global disruptions (including those in 1a) in the last 540 million years. The duration $$\tau$$ is the time over which $$m$$ grows. 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. Figure 1a adapted from [7], 1b adapted from [4].
Observational data sets the stage [4]. Figure 1a displays a geochemical time series containing two disruptions that occurred about 54 million years ago. The sharp downward pulses represent increases $$\Delta m$$ in the mass $$m$$ of inorganic carbon in the ocean. The ratio $$\Delta m/m$$ can be estimated from the amplitude and duration $$\tau$$ of the downswing. Figure 1b depicts results from 31 events over the last 540 million years. Roughly half of the events lie near the straight line. These characteristic events share a similar specific rate $$r_c$$ in the relation $$\Delta m/m= r_c \tau$$. Although $$r_c$$ is biogeochemically significant [4], it may also be dynamically significant because it appears to separate four of the five mass extinction events from nearly all other disruptions.

Carbon cycle disruptions are usually interpreted as proportionate responses to perturbations, such as enhanced emissions of volcanic CO2, extraordinary releases of methane, or changes in the rates at which organic carbon is sequestered in rock. But such a variety of stressors would seem unlikely to exhibit a common specific rate. Characteristic events may instead reflect the intrinsic dynamics of the carbon cycle.

A simple model of the marine carbon cycle illustrates how this could work [5]. We consider the upper ocean to be a well-mixed chemical reactor, open to an incoming flux $$j_i$$ of dissolved calcium carbonate from rivers and a maximum outgoing flux $$b j_i$$ to sediments. Dissolved inorganic carbon in the ocean takes several forms; it suffices to consider only two. We track the concentration of total dissolved inorganic carbon (carbonate and bicarbonate ions in addition to CO2), denoted by $$w$$, and the concentration of carbonate ions $$\textrm{CO}^{2-}_{3}$$), denoted by $$c$$. We also allow for external input of CO2 at rate $$\nu j_i$$.

One important feedback concerns the outgoing flux. Above a critical concentration $$c_p$$, carbonate minerals are preserved in sediments; below $$c_p$$, they dissolve. Averaged over the oceans, the transition is smooth rather than sharp. We therefore assume that the outflux is proportional to a sigmoidal function $$(s(c,c_p))$$ that grows from zero to one as $$c$$ increases to values well above $$c_p$$.

Another feedback concerns the carbonate system’s interaction with planktonic organisms. When CO2 is added to seawater, some of it combines with carbonate ions to form bicarbonate ions; this lowers the concentration $$c$$ of carbonate ions. If $$c$$ becomes sufficiently small—less than a concentration $$c_x$$—planktonic organisms with calcium carbonate shells fare poorly. Their shells normally provide “ballast” that causes them to sink to the deep sea with detrital organic carbon. But if fewer calcifying organisms exist in the upper ocean, less organic carbon sinks to the seafloor and more is respired to CO2 near the sea surface, causing upper ocean CO2 levels to increase even further. We express this  process by a constant $$\theta$$ times a sigmoidal function $$\bar{s} (c,c_x)$$, where $$\bar{s} = 1-s$$. Finally, we assume that the concentration of total dissolved inorganic carbon tends to relax toward a concentration $$w_0$$ at a timescale $$\tau_w$$. Today, $$\tau_w \sim 10^4$$ years. These considerations yield the carbon cycle model:

$\dot{c}/f(c) = j_i [ 1 - b s(c,c_p) - \theta \bar s(c,c_x) - \nu] + (w - w_0)/\tau_w \tag1$

$\dot{w} = j_i [ 1 - b s(c,c_p) + \theta \bar s(c,c_x) + \nu ] - (w - w_0)/\tau_w. \tag2$

The function $$f(c)$$ approximates the manner in which the carbonate system “buffers” the addition or removal of chemical species.

Depending on the parameters, the model exhibits a stable fixed point, a stable limit cycle, or both. Here we focus on the stable fixed point’s response to perturbation when the fixed point is near the bistable region. We imagine that the system is initially prepared with the CO2 injection parameter $$\nu=0$$, but we set $$\nu=\nu_0>0$$ for all times $$t \geq 0$$. Figures 2a and 2b show that the system spirals toward a new steady state. Yet when $$\nu_0$$ is larger than a threshold $$\nu_c$$, the system undergoes a large excitation before approaching the new fixed point (see Figures 2c and 2d). The model is therefore excitable, similar to models of action potentials in neurons [2]. The excitation corresponds to transient ocean acidification, since the addition of CO2 reduces not only $$c$$ but also the pH.

Figure 2. Phase plane and time series representations of perturbations to the carbon cycle model below (2a and 2b) and above (2c and 2d) the excitation threshold. Image courtesy of [5].

Because the size and timescale of excitations are properties of the system rather than its perturbation, an excitable carbon cycle could explain several features of the data in Figure 1b [5]. Characteristic events, which fall close to the straight line, would represent near-threshold excitations. The mass extinction events above the line would be associated with perturbations that significantly exceed the threshold. Events below the line may simply represent slow quasistatic change.

In the real carbon cycle, initiation of an excitation by CO2 injection would likely last for only a finite time $$t_i$$, e.g., $$\nu = \nu_0$$ for $$0\leq t \leq t_i$$. In the carbon cycle model, the excitation threshold $$\nu_c$$ is independent of $$t_i$$ when $$t_i$$ exceeds approximately one damping time $$\tau_w$$. But if $$t_i \ll \tau_w$$, perturbations can be damped before they lead to excitation. In this case, the threshold depends on the total mass injected and $$\nu_c \propto t_i^{-1}$$. Figures 3a and 3b illustrate how this works. Similar phenomena exist in other excitable systems [2, 8].

Figure 3. The excitation threshold $$\nu_c$$ depends on the duration of forcing. 3a. The response $$w(t)$$ to CO2 injection (in 3b) at varied dimensionless rates $$\nu$$ over different durations $$t_i$$. 3c. Comparison of perturbations of the modern and end-Cretaceous [6] carbon cycles to a hypothesis for the upper bound of the excitation threshold. Four projections of the modern perturbation are depicted [5]. The end-Cretaceous case $$\overline{KT}$$ represents the upper half of the confidence interval for CO2 emitted by massive volcanism tens of thousands of years before the end-Cretaceous extinction [6]. Because the excitation threshold scales like $$t^{-1}$$, both the modern and ancient perturbations are equivalently near it. Figure 3a and 3b courtesy of Daniel Rothman, 3c adapted from [5].

The flux of CO2 entering today’s oceans is nearly two orders of magnitude greater than the growth rate of characteristic events, which provides a rough estimate of the threshold’s upper bound under a step-like injection. Yet the centennial timescale of the modern perturbation is about two orders of magnitude smaller than today’s damping timescale. The shorter injection timescale cancels out the stronger injection rate, and the modern influx lies near the hypothesized upper bound rather than greatly exceeding it.

To understand what this might mean, consider the carbon cycle’s behavior before the end-Cretaceous extinction (and the dinosaurs’ demise). Although the extinction is widely attributed to a bolide impact, modern geochronological methods reveal a roughly 104-year pulse of massive volcanism tens of thousands of years prior to impact [6]. The resulting CO2 injection occurred at an estimated peak rate that is approximately one percent of the maximum projected 21st-century mean flux to the oceans [5]. Both the modern and end-Cretaceous CO2 fluxes therefore lie near the threshold’s upper bound (see Figure 3c). Massive volcanism associated with the end-Permian and end-Triassic extinction events may be similarly pulsed.

The upshot is twofold. First, modest perturbations of the carbon cycle beyond a threshold may have led to significant disruptions of the ancient Earth system, possibly including mass extinction. Second, today’s strong perturbation appears to be near an equivalent threshold. Are we therefore headed for a sixth extinction? The dynamical system of equations $$(1)$$ and $$(2)$$ is a toy model; among other limitations, its assumption of a well-mixed ocean neglects possible feedbacks on timescales less than approximately 103 years. Yet the hypothesis of excitability is reasonable and the $$t^{-1}$$ scaling law provides something novel and important: a means of rescaling the past to the present. Analysis of alternative models and acquisition of more data will help test these ideas. The risk of catastrophe makes fundamental progress imperative.

Acknowledgments: I thank Constantin Arnscheidt for helpful discussions. This work was supported by NASA astrobiology grant NNA13AA90A and NSF grant 1338810.

References
[1] Archer, D. (2010). The Global Carbon Cycle. Princeton, NJ: Princeton University Press.
[2] Izhikevich, E.M. (2007). Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. Cambridge, MA: MIT Press.
[3] Rothman, D.H. (2015). Earth’s carbon cycle: a mathematical perspective. Bull. Amer. Math. Soc., 52, 47-64.
[4] Rothman, D.H. (2017). Thresholds of catastrophe in the Earth system. Sci. Adv., 3(9), e1700906.
[5] Rothman, D.H. (2019). Characteristic disruptions of an excitable carbon cycle. PNAS, 116, 14813-14822.
[6] Schoene, B., Eddy, M.P., Samperton, K.M., Keller, C.B., Keller, G., Adatte, T., & Khadri, S.F.R. (2019). U-Pb constraints on pulsed eruption of the Deccan Traps across the end-Cretaceous mass extinction. Science, 363, 862-866.
[7] Westerhold, T., Röhl, U., Frederichs, T., Agnini, C., Raffi, I., Zachos, J.C., & Wilkens, R.H. (2017). Astronomical calibration of the Ypresian timescale: implications for seafloor spreading rates and the chaotic behavior of the solar system? Clim. Past, 13, 1129-1152.
[8] Wieczorek, S., Ashwin, P., Luke, C.M., & Cox, P.M. (2011). Excitability in ramped systems: the compost-bomb instability. Proc. Roy. Soc. A, 467(2129), 1243-1269.

Daniel H. Rothman is a professor of geophysics and co-director of the Lorenz Center in the Department of Earth, Atmospheric, and Planetary Sciences at the Massachusetts Institute of Technology.