discussion forum

Thresholds of catastrophe in the Earth system

Daniel H. Rothman, Science Advances, Sep 25 2017 - 08:45

The history of the Earth system is a story of change. Some changes are gradual and benign, but others, especially those associated with catastrophic mass extinction, are relatively abrupt and destructive. What sets one group apart from the other? Here, I hypothesize that perturbations of Earth’s carbon cycle lead to mass extinction if they exceed either a critical rate at long time scales or a critical size at short time scales. By analyzing 31 carbon isotopic events during the past 542 million years, I identify the critical rate with a limit imposed by mass conservation. Identification of the crossover time scale separating fast from slow events then yields the critical size. The modern critical size for the marine carbon cycle is roughly similar to the mass of carbon that human activities will likely have added to the oceans by the year 2100.


Five times in the Phanerozoic [the past 542 million years (My)], more than three-fourths of marine animal species have vanished in mass extinctions (1). Each of these events is associated with a significant change in Earth’s carbon cycle (2, 3). However, there are a great many disturbances of the carbon cycle (3, 4) whose imprint in the geologic record is almost exclusively environmental. These other events are relatively benign, with extinction rates barely, if at all rising, above background levels (57).

Two well-studied examples illustrate these distinctions. The end-Permian extinction [~252 million years ago (Ma)], the most severe mass extinction in the Phanerozoic (8), plays out over a period of 104 to 105 years; the extinction interval immediately follows a perturbation of the carbon cycle of similar duration (9). The Paleocene-Eocene Thermal Maximum (~55.5 Ma) is a carbon cycle event of roughly similar time scale, with unambiguous signs of global warming and ocean acidification (10, 11). It is associated with a significant extinction of benthic foraminifera (12, 13), but there is no evidence of other major extinctions (11, 13, 14). What makes these two events so different? The unambiguous co-occurrence of the end-Permian extinction with massive volcanism (15) provides one indication, as does the greater magnitude of the end-Permian carbon isotopic excursion compared to that of the Paleocene-Eocene Thermal Maximum (10, 16). However, no general distinction, applicable throughout the Phanerozoic, exists (57).

The idea that mass extinction results from catastrophic environmental change was largely developed more than two centuries ago by Cuvier (17). But what defines such an environmental catastrophe? Newell (18) implicitly provides a practical definition: “Extinction… is not simply a result of environmental change but is also a consequence of failure of the evolutionary process to keep pace with changing conditions in the physical and biological environment” (18). In other words, catastrophe results when the rate of environmental change exceeds a threshold.

To consider such ideas within the context of the carbon cycle, suppose that the mass m of inorganic carbon in the oceans is increased from its steady-state value, m*, by an amount Δm over a duration τenv of time. In terms of the relative change M = Δm/m*, Newell’s statement hypothesizes that mass extinction results when Menv exceeds a critical rate rc. However, one recognizes immediately that this condition cannot apply for impulsive change as τenv → 0. In this case, M must instead exceed a finite critical size Mc. Putting it all together, catastrophe follows whenEmbedded Imagewhere τx is the crossover time scale separating fast from slow change. Because the asymptotic regions 1 and 2 intersect at τenv = τx, the three unknowns Mc, rc, and τx can be assumed to satisfyEmbedded Image(3)In principle, Eq. 3 shows how the record of slow environmental change at geologic time scales informs our understanding of rapid change at human time scales, even though mechanisms at the different time scales may differ. Here, I provide an analysis of the geological record that makes contact with this picture.


I begin by constructing an empirical database of significant global perturbations of Earth’s carbon cycle during the Phanerozoic. Each perturbation is considered to be an isolated event, or excursion, in the evolution of the isotopic composition δ1 of inorganic (carbonate) carbon. Except for the most recent event (beginning at the Last Glacial Maximum), the analysis focuses exclusively on isotopic changes that decrease δ1. Each isotopic event is parameterized as shown in Fig. 1A.

Fig. 1 Parameterization and temporal distribution of carbon isotopic events in the database.

(A) Schematic diagram illustrating how the parameters δ1(0) and δ1env) are derived from an isotopic time series δ1(t). The time scale τenv represents the duration of time during which the isotopic composition δ1 decreases. (B) Histogram depicting the distribution of events in time.

The database contains each event’s geologic age, the time τenv over which δ1 decreases, and the isotopic compositions δ1(0) and δ1env) that bracket the change from time t = 0 to t = τenv. Whenever possible, estimates of τenv are obtained from isotope geochronology or astrochronology. The compilation contains a total of 31 events distributed throughout the Phanerozoic (Fig. 1B), including each of the so-called Big Five mass extinctions. The resulting data are contained in Table 1. Further information is available in Materials and Methods.

Table 1 Event names, event abbreviations, and data used to construct Figs. 2 and 3.

Values of δ1 and Δr are given in per mil (‰) and derive from the formula δ1 = 1000 × (RRstd)/Rstd, where abundance ratios R = 13C/12C are obtained from samples of inorganic (carbonate) carbon and Rstd is a standard ratio. The variables Δr, τenv, φ, and M are defined in the text. In Fig. 2, the six Eocene and four Miocene events are denoted by Eo and Mio, respectively.

View this table:

Figure 2A displays the isotopic response Δr = δ1env) − δ1(0) of each event as a function of its time scale τenv. Two features of the data immediately stand out. First, there is a relative absence of large excursions at short time scales, suggesting that there may be a characteristic limitation to the rate at which the marine inorganic carbon reservoir can change its isotopic composition. Second, four of the Big Five mass extinctions appear to exhibit anomalously fast rates of change. The fifth mass extinction, in the late Devonian at the Frasnian-Famennian boundary, is clearly different, lying in the slow region in the lower right-hand corner (labeled “FF”). McGhee (19) and Bambach et al. (20) have suggested that the Frasnian-Famennian event is not an extinction event but instead represents a decrease in the rate at which new species originate. Its relatively small and slow isotopic perturbation provides indirect support for this idea.

Fig. 2 Depictions of the size and time scale of carbon isotopic events.

(A) Magnitude of the isotopic shift Δr = δ1env) − δ1(0) as a function of its duration of time, τenv. (B) Dimensionless mass perturbation M = |Δm|/m* as a function of the dimensionless time scale φcτenv/2τ0 for each of the events depicted in (A). The straight (identity) line denotes the equality predicted by Eq. 5 when φ = φc = 0.23 ± 0.07. Event abbreviations are defined in Table 1. Error bars (Materials and Methods) are guides for interpretation.

To understand these observations more deeply, I transform the geochemical changes to physical fluxes (Materials and Methods). The transformation assumes that the global carbon cycle is composed of two reservoirs: inorganic and organic carbon. Carbon is exchanged between the two reservoirs by photosynthesis and respiration. Inorganic carbon of isotopic composition δi enters the cycle with flux ji. Both inorganic and organic carbon exit the cycle when buried as rock.

The downturn of δ1 is assumed to be driven by a time-dependent flux Embedded Image of inorganic carbon with isotopic composition Embedded Image, which is set equal to the average Phanerozoic isotopic composition of organic carbon [− 27.8‰ (21)]. In other words, Embedded Image represents an imbalance in the loop between photosynthesis and respiration on time scales up to τenv. Its origin is unspecified, but its growth is assumed to be linear. Taking each event to begin at t = 0 and defining the normalized flux Embedded Image, where Embedded Image is the unperturbed input flux, one then hasEmbedded Image(4)where φ is an unknown, event-specific, dimensionless constant that specifies the maximum normalized flux Jenv). An alternative model could be a constant perturbation J = φ. But such a choice would result in equilibration to a new steady state in the many instances in which τenv is greater than the turnover time Embedded Image (thousand years) (Materials and Methods). The analysis here instead assumes that perturbations grow until t = τenv. The linear model (Eq. 4) can therefore be viewed as a first-order approximation of a flux function that satisfies this condition. The normalized mass M in Eqs. 1 and 2 derives from its integral:Embedded Image(5)

To determine the flux J(t) consistent with each isotopic event, I calculate the value of φ that solves differential equations for the evolution of δ1(t), M(t), and J(t) given the initial values Embedded Image, J(0) = 0, and M(0) = 0 and the final value Embedded Image. The values of Embedded Image, Δr, and τenv are obtained for each event from the database. The solution of this boundary value problem for φ provides the evolution of δ1(t) that matches its observed initial and final states while remaining consistent with differential equations describing changes in the isotopic composition of marine inorganic carbon in response to the flux J(t). It also provides, via Eqs. 4 and 5, the rate at which J grows and thus the mass M.

Figure 2B shows M as a function of τenv, rescaled according to Eq. 5, for each event. The gross features seen in the raw data of Fig. 2A remain, but now, the appearance of a limiting rate, suggested once again by the relatively barren upper-left half of the plot, is sharper. Moreover, the limiting rate now has a physical interpretation. From Eq. 4, one sees that the normalized flux perturbation reaches its maximum value, J = φ, when t = τenv. The limiting rate therefore represents a critical dimensionless flux φc beyond which the carbon cycle rarely ventures. When it does, mass extinction often follows.

I hypothesize that φc derives from organic carbon that would normally be buried as rock but is instead converted to CO2 and remobilized into the marine carbon cycle. The rate of this remobilization can be no greater than the rate at which organic carbon is immobilized (that is, buried) in an unperturbed steady state. During the Phanerozoic, this limiting rate equals 23 ± 7% of the steady input (or output) flux Embedded Image (Materials and Methods). Therefore, φc = 0.23 is the maximum possible normalized flux perturbation in a mass-conserving carbon cycle with no anomalous sources or sinks. The straight line in Fig. 2B, which shows Mc, τenv) as defined by Eq. 5, represents this hypothesis quantitatively.

This line separates two kinds of carbon cycles. Below the line, perturbations can be explained by a partial cessation of organic carbon burial. Above the line, they cannot. On the line, the small leak represented by burial is temporarily closed when the perturbation reaches its maximum strength at t = τenv. Roughly half of the events in the database fall within the error bars of this boundary. Four of the Big Five mass extinctions lie above the upper error bar, whereas nearly all other events fall below it. These observations suggest that noncatastrophic perturbations represent flux changes along a fixed set of pathways, whereas catastrophic perturbations result from the introduction of new pathways acting as anomalous sources or sinks. I refer to the borderline perturbations as critical events. The preponderance of critical events supports the idea that large yet benign flux perturbations tend to grow until their source—the small fraction of organic matter that would normally be destined for burial—is fully depleted via respiration.

The critical events are likely driven to their limit by the intrinsic dynamics of the carbon cycle over a time scale set by the characteristic relaxation time associated with a negative feedback. However, the variation of τenv over 104 to 106 years would seem to rule out a unique damping mechanism. This disqualification would necessarily hold if the carbon cycle did not change over geologic time. But evidently it did: Fig. 3A shows that the upper envelope of event size M decreases with time. To understand this evolution, note that Eq. 5 predicts that changes in the size of critical events, for which φ = φc, derive directly from changes in τenv. Figure 3B shows that this time scale—which should be roughly equivalent to the damping time—decreases toward the present.

Fig. 3 Evolution of size and time scale.

(A) Dimensionless mass perturbation M = |Δm|/m* as a function of the geologic age of the event. The relative absence of smaller events in the deeper past likely derives from poor data and a lack of interest. The absence of large events in the more recent past, however, cannot be explained by this bias. (B) Time scale τenv of critical events (events falling within the error bars of the straight line in Fig. 2B) as a function of geologic age. Error bars (Materials and Methods) are guides for interpretation.

The time scale of critical events appears to have declined since the early Triassic (~250 Ma) by about one order of magnitude, to about 20 to 30 ky since the Miocene (~16 Ma). The smaller time scale is similar to the ~10-ky time scale of the modern oceans’ homeostatic response to a change in pH (2225). The time scale’s decrease since the Triassic is consistent with the evolutionary expansion of the modern biological pump beginning approximately 220 Ma (26, 27), corresponding to the earliest observations of coccolithophores and associated calcareous phytoplankton in the fossil record (28, 29). As the biological pump strengthened, the export of precipitated calcium carbonate from surface waters to deep ocean sediments increased, leading to growth of the deep-sea sedimentary carbonate reservoir (26, 27, 30). When pH decreases in modern oceans, dissolution of sedimentary carbonate minerals is one of the principal processes via which these disturbances are damped (31, 32). Because this mechanism strengthened over the last 220 My, critical perturbations of the marine carbon cycle would have become shorter in time and smaller in magnitude, but not necessarily less dangerous, because their maximum flux would have remained φc. These observations are each consistent with the suggestion that the carbon cycle became less variable after the introduction of the deep-sea carbonate sink (26, 27, 30, 33). They also identify the crossover time scale τx with the homeostat time scale, implying τx ~ 10 ky (2225) in modern oceans.

With these results in hand, the modern critical size Mc of an asymptotically fast perturbation of the marine carbon cycle now follows. First, insertion of Mc, τenv) from Eq. 5 into Eq. 2 provides the critical rate rc = φc/2τ0. Equation 3 then predictsEmbedded Image(6)Embedded Image(7)where, as a guide, a total uncertainty of ±50% is assumed. The modern oceans contain about 38,000 Pg C; the critical dimensional mass is therefore about 310 ± 155 Pg C.

From 1850 to the present, human activities have resulted in the addition of about 155 ± 25 Pg C to the oceans (34). Projections for further carbon uptake depend strongly on the trajectory of fossil fuel emissions and land use, among other factors (34, 35). Figure 4 compares the critical mass to the present accumulation and four projections to 2100 obtained from coupled climate–carbon cycle models (34). The strictest emission scenario results in oceanic carbon uptake whose mean falls just below the critical mass; at the other extreme, the mean model uptake is about 74% greater than critical. Although the uncertainty of each prediction in Fig. 4 is considerable, all scenarios for cumulative uptake at the century’s end either exceed or are commensurate with the threshold for catastrophic change.

Fig. 4 Cumulative modern ocean uptake of carbon since 1850, up to the present (green) and projected to 2100 (blue), compared to the predicted critical mass of 310 Pg C (solid red line) and an assumed uncertainty of ±50% (dashed red lines).

Projections (34) are given for four representative concentration pathway scenarios RCPx, where x represents the radiative forcing, in units of W/m2, deriving from accumulated emissions in the year 2100 (35). At the extreme ends of the projections, RCP2.6 represents the range of lowest-emission scenarios in the scientific literature, and RCP8.5 represents the high range of nonclimate policy scenarios. Of the two intermediate pathways, RCP4.5 corresponds to an emission pathway resulting from many climate policies found in the literature, whereas RCP6.0 is representative of most scenarios without limitations on emissions (35). The present cumulative uptake is obtained by adding 6 years of an annual uptake rate of 2.3 Pg C year−1 (34) to the 2011 total of Ciais et al. (34).


The value of the empirical results in Figs. 2 and 3 lies not in any individual data points but rather in their collective trends. The following picture emerges. Carbon isotopic events rarely exceed a maximum isotopic shift that grows roughly like the logarithm of their time scale. This upper bound appears related to the minimum rate—zero—at which organic carbon can be immobilized as rock. Events outside this limit result from a fundamental disturbance of the carbon cycle, possibly related to unstable dynamics, mass extinction, or both.

These conclusions follow from analyzing all isotopic events the same way. Exceptions are, however, expected. For example, four events (Ediacaran-Cambrian, Nemakit-Daldynian-Tommotian, Paleocene-Eocene Thermal Maximum, and Miocene Climatic Optimum 2) unaccompanied by mass extinction exceed the upper error bar of the critical rate. If, say, these events were driven by dissociation of methane hydrates [for example, (36)] rather than respired organic carbon, the isotopic composition of the perturbative flux would be much lighter (about −60‰). The estimated size of these events would then drop by more than 50%, and they would each lie below the critical line. This observation may help explain the Paleocene-Eocene Thermal Maximum’s modest biotic impact. The Frasnian-Famennian extinction provides another exception. Supposing that it is indeed a mass extinction, its presence well below the critical line illustrates an important point: Mass extinctions need not be caused by disruptions of the carbon cycle (2).

Modern investigations of mass extinctions often emphasize a plurality of causes. Erwin’s “complex web of causality” (8, 37) addresses how a combination of volcanism, climate change, marine anoxia, methane release, and other environmental stressors may have contributed to the end-Permian extinction. Recent studies of the end-Cretaceous extinction consider massive volcanism (38) in addition to a bolide impact (39). Flood basalt eruptions are also clearly associated with the end-Triassic (40) and end-Permian (15) extinctions, but their contribution to CO2 levels is ostensibly modest (41). Evidently, the carbon cycle both indicates and excites Earth system change. These dual roles merge, however, if external perturbations cause the cycle to respond by magnifying the initial disturbance. System-wide instability may then follow. Because the critical rate rc bounds qualitatively different dynamical regimes, perturbations that exceed rc (at time scales much greater than τx) suggest such unstable evolution. The carbon cycle thus becomes one of many environmental stressors, and an array of causes is naturally implicated.

Multiple causes may also be responsible for a “sixth extinction” (1, 42). However, the anthropogenic disturbance of the carbon cycle merits its own appraisal. If an anomalous flux greater than Mcx brings the carbon cycle past a stability threshold at time scales much greater than τx ~ 104 years, then the instantaneous addition of an anomalous mass greater than Mc will excite similar behavior. The existence of these thresholds is this paper’s central hypothesis. If they exist, one would expect that an instability resulting from an impulsive addition of CO2 would play out over a time scale with the same order of magnitude as the damping time scale τx. The upshot is that an unstable trajectory would reach its maximum extent roughly 104 years after the threshold is breached. But how that process plays out remains unknown.


The event database

After a brief discussion of possible sources of error, this section summarizes how each event has been identified and quantified. Geological, paleoenvironmental, and paleobiological context is also given, as are relevant references. All data, including the identification of abbreviations, are provided in Table 1.

When multiple observations are used to compute a mean, the sample SD, the unaveraged observations, or both are noted in the event’s description. However, these quantities only indicate the variability of a particular set of observations rather than an estimate of the observational error. As with any compilation of isotopic data [for example, Hayes et al. (43)], inadequate sampling may be a more significant problem. Here, it takes two forms: the sheer paucity of observations, and biases deriving from not only the original investigators—who may favor larger isotopic excursions—but also the needs of the present compilation, which favors carbon isotopic studies associated with strong geochronological constraints. An additional source of error derives from incompleteness of the sedimentary record. For example, a hiatus in sediment deposition, or erosion after deposition, may make the true carbon isotopic minimum effectively unavailable to the modern observer, thereby leading to an underestimate of |Δr|.

The ideal geochronological control would provide absolute dates at the initiation and termination of isotopic excursions. Only one event in the database—at the Ediacaran-Cambrian boundary—meets this standard. Among the others, all but 1 of the 18 events extending from the Toarcian (182 Ma) to the Miocene (16 Ma) are timed by astrochronology, including, in some cases, information obtained from isotope geochronology. The duration of five of the older, pre-Jurassic events is determined from at least two dates provided by U-Pb geochronology that bracket less than twice the accumulation of sediment deriving from the event itself; a sixth event (end-Permian) comes close to this standard.

The remaining pre-Quaternary events are timed by biostratigraphy with reference to the Geologic Time Scale (44); their time scales are therefore the most poorly constrained. Two of them, the end-Ordovician and Frasnian-Famennian events, are among the Big Five mass extinctions. A third mass extinction—the end-Triassic—is constrained by both isotope geochronology and astrochronology, but these constraints are not directly associated with the carbonate–carbon isotope excursion. If these three mass extinction events and the other four events with poor temporal constraints (Cambrian Spice, Tournaisian, mid-Capitanian, and Albian-Cenomanian) are removed from the database, the trends in Figs. 2 and 3 remain intact and no new trends appear.

Nevertheless, the unknown errors in the observations render quantitative error analysis infeasible. The representative error bars for Δr, M, and τenv in Figs. 2 and 3 are instead guides for interpretation.

The following list proceeds from the oldest to the youngest event.

Ediacaran-Cambrian boundary, ~541 Ma. The magnitude and duration of this negative isotopic excursion derive from the Oman carbon isotope data of Amthor et al. (45) and the U-Pb geochronology of Bowring et al. (46). Other data from Morocco, Siberia, Mongolia, and China (47, 48) suggest global consistency. The absence of the Ediacaran biota in the Cambrian has led to the suggestion that the Ediacaran biota vanished by mass extinction at the end of the Ediacaran, possibly related to the environmental event signaled by this excursion (45, 4952). However, gradual extinction coupled with permanent environmental changes unfavorable to preservation of the Ediacaran biota may also be possible (53).

Nemakit-Daldynian-Tommotian boundary, ~525 Ma. Both the carbon isotopic excursion and the U-Pb dates establishing the duration of this Early Cambrian event are from the Moroccan data of Maloof et al. (47, 54); data from Siberia, Mongolia, and China confirm its global nature (47, 48). No mass extinction has been associated with this interval. However, diversity (number of genera) and disparity (number of classes) sharply increase in the Tommotian, primarily because of the rise of small shelly fossils (55). The Nemakit-Daldynian-Tommotian boundary is also associated with a transition from seawater chemistry that favored aragonite precipitation to seawater favoring calcite precipitation (47, 56).

Cambrian Spice, ~497 Ma. The duration and magnitude of this late Cambrian event, also known as the Steptoean positive carbon isotope excursion, are taken from the Australian data of Saltzman et al. (57). The event has also been found in North America, China, Kazakhstan, and Sweden (58, 59). The initiation of this positive isotopic excursion coincides with an extinction of trilobites (58).

End-Ordovician event, ~446 Ma. The latest stage of the Ordovician—the Hirnantian—simultaneously exhibits a positive isotopic excursion and one of the Big Five mass extinctions. I estimated the magnitude of the negative limb of the excursion from data obtained at Vinni Creek and Monitor Range, Nevada (60), Copenhagen Canyon, Nevada (61), Blackstone Range, Yukon (60), Wangjiawan, south China (62), the Kardla drill core, Estonia (63), and Anticosti Island, Quebec (64). The seven data sets yield a mean excursion of 5.6 ± 0.9‰. Each data set except for Copenhagen Canyon has been correlated to graptolite zones by Gorjan et al. (62). By estimating the fraction of the persculptus and extraordinarius zones in which the excursion occurs at each site and estimating those time intervals from Gradstein et al. (44), I found that the negative limb of the excursion lasted approximately 230 ky.

Silurian Mulde event, ~428 Ma. A sequence of two positive excursions, known as the “Mulde event,” occurs in the Homerian stage of the Silurian, during a period of graptolite extinctions known as the “Big Crisis” (65). The event is apparently global (65). Cramer et al. (66) have recently bracketed the earlier of the two excursions between two U-Pb dates obtained in Gotland, Sweden. I focused on the negative downswing, which occupies somewhat more than half of the dated interval, and estimated its duration to be about 260 ky. As shown by Cramer et al. (66), the Gotland excursion correlates well with observations in the West Midlands, England, with the West Midlands excursion having a magnitude of about 4‰ and the Gotland excursion about 2‰, a range that is consistent with other observations (65, 66). I therefore estimated the magnitude to be 3 ± 1‰.

Silurian Lau event, ~423 Ma. A significant global positive carbon isotopic excursion known as the “Lau Event” occurs in the Ludlow Epoch of the Silurian (65, 6772). The excursion typically peaks at 6 to 8‰ [for example, Lehnert et al. (71)], though values as small as 3 to 4‰ and as high as 10 to 11‰ have been observed in North America (68) and Sweden (70), respectively. Given that uncertainty, I adopted a conservative mean peak value of 7‰. The negative downswing after the peak typically bottoms out at about 1‰, resulting in a typical total negative shift of about 6‰. Cramer et al. (72) have recently bracketed the duration of the excursion by U-Pb geochronology, finding that it must be less than about 1.17 My. Following their analysis, I took the upswing and downswing to be roughly symmetric, each occupying about one-half of an approximately million-year interval, implying that the duration of the downswing is about 500 ky.

Frasnian-Famennian boundary, ~372 Ma. The boundary between the Frasnian and Famennian stages in the Late Devonian is associated with one of the Big Five mass extinctions (5, 6). But the drop in diversity at the boundary is thought to be a consequence of a lack of originations rather than an elevated extinction rate (19, 20). The boundary is associated with a global positive carbon isotopic excursion, known as the Upper Kellwasser Event (73). The magnitude and duration of the excursion are taken from the compilation of European Devonian data of Buggisch and Joachimski (74).

Tournaisian event, ~352 Ma. Buggisch et al. (75) presented carbon isotopic data from Belgium, Ireland, France, western Canada, and the western United States that collectively exhibit a pronounced positive isotopic excursion in the Tournaisian stage of the Carboniferous. Using their biostratigraphic time scale, I estimated the duration of the negative limb of the excursion to be 1.51 My; the magnitude of the excursion in their compilation is about 3.5‰.

Mid-Capitanian event, ~262 Ma. A mass extinction in the Capitanian Stage of the Middle Permian is associated with a negative isotopic excursion of carbonate carbon of about 5‰ over about 500 ky (76, 77). Observations of the event in south China suggest that the extinction event coincides with the initiation of Emeishan volcanism (78) and the onset of the excursion (76).

End-Permian extinction, ~252 Ma. The end-Permian extinction is considered the most severe of the Big Five (8). A significant negative isotopic excursion occurs at the time of the extinction, just below the Permian-Triassic boundary. Korte and Kozur (16) reviewed observations of the excursion in 40 localities worldwide and concluded that the excursion ranged from 4 to 7‰; I therefore took the magnitude of the excursion to be at the center of that distribution (5.5‰). The 60-ky duration of the excursion derives from the well-studied section at Meishan, south China, where Cao et al. (79) have provided high-resolution carbon isotopic data and U-Pb geochronology provides strong constraints on the time scale (9).

Early Triassic, ~251 Ma. Several significant isotopic excursions of unknown origin follow the end-Permian extinction (8083). Galfetti et al. (83) provided carbon isotopic data tied to U-Pb geochronology from the Loulou Formation, northwest Guanxi, south China. The only temporally constrained negative isotopic shift occurs in a negative excursion that begins approximately 251.22 Ma (83) in the mid-Smithian and declines approximately 3.5‰ over approximately 250 ky.

Triassic-Jurassic boundary, ~201 Ma. One of the Big Five extinctions, the end-Triassic event is temporally associated with the emplacement of the Central Atlantic magmatic province (40). Geochemically, the most well-studied section is at St. Audries Bay in southwest England, where the extinction interval coincides with a rapid fall and subsequent rise of the isotopic composition of organic matter (84) that lasts 20 to 40 ky according to astrochronology (85). This so-called initial excursion in organic carbon is widely observed (86). However, its significance for understanding the evolution of marine dissolved inorganic carbon is unclear (86, 87), partly because it may be associated with observed changes in flora (87)—and therefore variations in the organic matter itself—and partly because inorganic isotopic compositions need not track organic isotopic compositions (43). There are unfortunately few well-resolved isotopic studies of carbonate carbon associated with the end-Triassic event (8893). These carbonate studies also reveal a negative excursion associated with the extinction, but it is often unclear if they represent the “initial excursion” or the later 120-ky-long “main excursion” seen in organic carbon (84, 85). A recent review (94) cautions that diagenetic alteration may have corrupted carbonate data, yet a careful study of lithological effects at an Italian section suggests that the carbonate excursion may indeed be primary (87). An important additional problem is the need to estimate the time scale of the carbonate excursion. The data of Clémence et al. (92) provide a solution to these problems. In their analysis of inorganic and organic carbonate at the Tiefengraben section in the Austrian Alps, they found an initial 2‰ negative excursion in carbonate that extended to the second significant minimum in their organic data (92), which in turn correlates reasonably well with the astrochronologically calibrated organic data at St. Audries Bay (85). Their initial carbonate excursion was then found to last approximately 50 ky [which included one precessional cycle of the “pre-recovery” interval identified by Ruhl et al. (85)]. This 2.0‰ shift is also consistent with isotopic analyses of carbonate in pristine oysters located near St. Audries Bay (91).

Toarcian oceanic anoxic event, ~183 Ma. The early Toarcian oceanic anoxic event (95) is widely observed in Europe as a negative excursion in the isotopic composition of carbonate (9699). To estimate its magnitude and duration, I used the high-resolution carbonate record of Hesselbo et al. (96) obtained at Peniche (Portugal). The astrochronology of Suan et al. (97) found a duration of 150 ky (their interval C1) for the negative isotopic excursion, which is consistent with the subsequent analysis of Huang and Hesselbo (100). Using this chronology, recent U-Pb dating by Burgess et al. (101) found that the intrusive magmatism associated with the Ferrar large igneous province is synchronous with the negative excursion.

Aptian oceanic anoxic event, ~120 Ma. The negative isotopic excursion associated with the Aptian oceanic anoxic event has been observed worldwide in open ocean records of carbonate carbon (102105). My estimate of the magnitude and duration of the excursion derived from the astrochronology of Malinverno et al. (106) performed on the carbonate record of Erba et al. (103). The negative excursion spans the interval C3 [originally identified by Menegatti et al. (102)] and lasts approximately 47 ky.

Albian-Cenomanian boundary, ~100 Ma. This complex carbon isotope event occurred near the Albian-Cenomanian boundary (mid-Cretaceous) in several European sections (107). Here, I focused on the carbon isotopic data obtained at Ocean Drilling Program (ODP) Site 1050 at Blake Nose, western North Atlantic, as presented by Ando et al. (108). When plotted using the time scale of Huber et al. (109), the ODP data show an unambiguous 0.7‰ negative excursion across the boundary over approximately 110 ky.

Mid-Cenomanian event, ~95.9 Ma. The mid-Cenomanian positive carbon isotopic excursion has been clearly observed in European sections (107), the North Atlantic (110), western North America (111), and elsewhere. My analysis follows from the correlation of the “English Chalk” reference carbonate curve (107) to carbon isotopic data from the Cretaceous Western Interior Seaway of North America (111, 112). I used the time scale of Eldrett et al. (112), which is based on a combination of astrochronology and U-Pb geochronology. As indicated in figure 11 of Eldrett et al. (112), the negative swing of the excursion begins at about 96.5 Ma, followed by a decline of about 0.85‰ over roughly 138 ky.

Cenomanian-Turonian boundary event, ~94.2 Ma. Otherwise known as oceanic anoxic event 2, the positive isotopic excursion at the Cenomanian-Turonian boundary is probably the most widely observed of the Cretaceous isotopic events (107, 111, 112). As for the mid-Cenomanian event, my estimate of the time scale and magnitude of the Cenomanian-Turonian event follows from the correlation of the English Chalk reference carbonate curve (107) to carbon isotopic data from the Cretaceous Western Interior Seaway of North America (111, 112). I took the initiation of the event to be at the positive peak labeled B in figure 11 of Eldrett et al. (112) and found that the ensuing negative excursion declines by 1.60‰ in roughly 553 ky.

Cretaceous-Paleogene extinction, ~65.5 Ma. The end-Cretaceous negative isotopic excursion is associated with one of the Big Five extinctions (5), widely known for the extinction of dinosaurs and the Alvarez impact hypothesis (39). Its temporal association with the eruption of the massive Deccan volcanic province in India is less widely known (38). I obtained the magnitude and time scale of the carbon isotopic event from the deep-sea bulk carbon isotopic data obtained at ODP Sites 1210 (Northwest Pacific) and 1262 (South Atlantic), as presented by Alegret et al. (113) using an orbitally tuned time scale (114). I took the initiation and termination of the approximately 1.15 ± 0.03‰ (the mean of 1.16 and 1.13‰) negative isotopic excursion at each site to be where values of δ1 begin to change by at least 0.1‰, resulting in a mean event duration of roughly 26 ky (the average of 8.5 and 44 ky).

Early late Paleocene event, ~58.9 Ma. Also known as the mid-Paleocene biotic event, this negative excursion is synchronous with dissolution of carbonate and changes in the organization of benthic and planktonic ecosystems (115, 116). The magnitude and time scale of this negative excursion are taken from the high-resolution bulk isotopic record of Littler et al. (115).

Paleocene-Eocene Thermal Maximum, ~55.5 Ma. The negative isotopic excursion of the Paleocene-Eocene Thermal Maximum is perhaps the most studied carbon isotopic event in Earth history (11), in large part because of its association with significant climate warming and clear evidence of ocean acidification (10). The event is also associated with a significant extinction of benthic foraminifera (12); however, other groups of benthic and planktic microfossils show little or no extinction (11). The time scales of the isotopic excursion of bulk carbonate in two deep-sea cores, from ODP Site 1266 on the Walvis Ridge in the South Atlantic and ODP Site 690 in the Weddell Sea, Southern Ocean, have each been estimated independently by identification of orbital cycles (117) and the estimation of sedimentation rates from the concentration of 3He (118, 119). The mean of the resulting four estimates [corresponding to the cumulative time between the onset of the excursion and the termination of its “core” as summarized in table 1 of Murphy et al. (119)] is 83 ± 23 ky. I obtained the magnitude of the isotopic excursion, 2.7 ± 1.1‰, from the mean of 33 published analyses of bulk Paleocene-Eocene Thermal Maximum carbonates reviewed by McInerney and Wing (11). The initial value of the excursion was obtained similarly (11). Although these averages lack prejudice, they may nevertheless underestimate the excursion’s size (120) and overestimate its time scale (121).

Eocene Thermal Maximum 2, ~53.7 Ma. The Eocene is punctuated by several “hyperthermal events,” each represented by a negative excursion in the isotopic composition of carbonate carbon and dissolution of deep-sea marine carbonates. Eocene Thermal Maximum 2, one such event, follows the Paleocene-Eocene Thermal Maximum (also known as Eocene Thermal Maximum 1) by roughly 2 My. My estimate of the magnitude of Eocene Thermal Maximum 2 derived from averaging estimates from four high-resolution bulk isotopic records presented by Stap et al. (122). The records derive from ODP Sites 1262, 1263, 1265, and 1267, corresponding to paleowater depths ranging from about 1500 to 3600 m (122); the shallowest site yields an excursion of about 1.5‰, whereas the others are each about 1.0‰ (122), yielding a mean of 1.13 ± 0.25‰. The time scale is about 45 ky.

Eocene Hyperthermal H2, ~53.6 Ma. The magnitude and time scale of this Eocene hyperthermal event come from ODP Sites 1263, 1265, and 1267. Following Stap et al. (122), I took the excursion of bulk carbon isotopes to be about 0.6‰ and the time scale to be about 33 ky.

Eocene Hyperthermal I1, ~53.2 Ma. I obtained the magnitude (0.72‰) and time scale (41 ky) of this Eocene hyperthermal event from the high-resolution bulk carbon isotopic record of Littler et al. (115), which derives from ODP Site 1262.

Eocene Hyperthermal I2, ~53.1 Ma. This Eocene hyperthermal follows Eocene Hyperthermal I1 by about 100 ky. Using the same source (115) as for Eocene Hyperthermal I1, I found an excursion of about 0.61‰ and a time scale of 40 ky.

Eocene Thermal Maximum 3, ~52.5 Ma. I estimated the magnitude and time scale of this hyperthermal event from the high-resolution benthic carbon isotopic record of Lauretano et al. (123), obtained at ODP Sites 1262 and 1263. The two records are similar, yielding an excursion of about 0.8‰ and a time scale of about 37 ky.

Early Oligocene Event, ~33.5 Ma. The positive carbon isotopic excursion just above the Eocene-Oligocene boundary is associated with the initiation of permanent Cenozoic ice sheets on Antarctica (124126). I focused on the negative downswing to lighter isotopic compositions that followed. The time scale and isotopic change are derived from the high-resolution data collected at ODP Site 1218, using the astrochronological time scale presented by Coxall et al. (125). The negative shift extends over roughly 0.6 ± 0.1‰ during a period of about 430 ± 100 ky. The three ODP records presented by Zachos and Kump (126) are consistent with these estimates.

Miocene Climatic Optimum 1, ~16.9 Ma. Holbourn et al. (127) provided high-resolution carbon isotopic records spanning most of the Miocene Climatic Optimum, a period of warm climates ranging from 17.0 to 14.7 Ma that interrupted the longer-term trend of Cenozoic cooling. The records are obtained from Integrated ODP Site U1337 in the eastern equatorial Pacific Ocean, at a paleowater depth of 3500 to 4000 m. The Miocene Climatic Optimum began with a negative isotopic excursion, which I designated Miocene Climatic Optimum 1, that has also been identified in lower-resolution records from three other sites, elsewhere in the Pacific and in the Southern Ocean (127). I obtained the magnitude (0.5‰) and time scale (28 ky) from the bulk carbonate record. The event is synchronous with shoaling of the carbonate compensation depth (CCD).

Miocene Climatic Optimum 2, ~16.4 Ma. Holbourn et al. (127) identified three other Miocene Climatic Optimum negative excursions, which I designated Miocene Climatic Optimum 2 to Miocene Climatic Optimum 4, during which δ1 decreases sharply and the CCD appears to shoal. Using the same high-resolution record as for Miocene Climatic Optimum 1, I found that Miocene Climatic Optimum 2 has a time scale of 22 ky and a magnitude of 0.76‰.

Miocene Climatic Optimum 3, ~16.0 Ma. See Miocene Climatic Optimum 2. For Miocene Climatic Optimum 3, I found a time scale of 32 ky and a magnitude of 0.66‰. This event also appears in the records of ODP Site 1146 in the western Pacific Ocean (128) and ODP Site U1338 in the eastern equatorial Pacific (129).

Miocene Climatic Optimum 4, ~15.6 Ma. See Miocene Climatic Optimum 2. For Miocene Climatic Optimum 4, the excursion has a magnitude of 0.6‰ and a time scale of 29 ky.

Last Glacial Maximum–to–Holocene deglacial, ~0.021 Ma. Using a compilation of 480 records of benthic foraminiferal carbon isotope analyses, Peterson et al. (130) estimated a whole-ocean increase in the isotopic composition of dissolved inorganic carbon of 0.34 ± 0.19‰ from the Last Glacial Maximum to the Holocene. I took the time scale to be 18 ky, given their Last Glacial Maximum averages over 19 to 23 ky and Holocene averages over 0 to 6 ky. This is the only positive isotopic shift in the database in which the magnitude and time scale derive from the upswing; for the others, such as the Ordovician, the calculations correspond to the ensuing downshift.

The boundary value problem

The boundary value problem derives from a model of the global carbon isotopic system presented by Rothman et al. (131, 132). This model partitions carbon into two reservoirs. The inorganic reservoir contains a mass m of carbon with isotopic composition δ1; the organic reservoir contains carbon of isotopic composition δ2. In the Phanerozoic, the isotopic composition of the organic reservoir typically equilibrates with respect to inputs and outputs at time scales much faster than the inorganic reservoir (131). The evolution of δ1 is then described by equation 8 of Rothman et al. (132):Embedded Image(8)where ji is the flux of carbon of isotopic composition δi into the surface carbon cycle (deriving typically from volcanism and continental weathering), b2 is the burial flux of organic carbon out of the surface cycle, and ε = δ1 − δ2 is the isotopic fractionation between inorganic and organic carbon.

Suppose now that a new source of carbon with isotopic composition Embedded Image flows into the carbon cycle with flux Embedded Image, and that the burial flux b2 increases by an amount Embedded Image while ε remains constant. Equation 8 then becomesEmbedded Image(9)where Embedded Image, Embedded Image, and Embedded Image represent unperturbed quantities that satisfy the steady-state relationEmbedded Image(10)where Embedded Image is the steady-state isotopic composition of the inorganic reservoir. Subtraction of Eq. 10 from Eq. 9 then yieldsEmbedded Image(11)Here, Embedded Image is set to the average isotopic composition of Phanerozoic organic carbon (−27.8‰ (21)]. Consequently, the terms in Eq. 11 proportional to Embedded Image and Embedded Image are practically indistinguishable so that a positive perturbation of the input flux is equivalent to a negative perturbation of the burial flux. I therefore set Embedded Image and reinterpreted Embedded Image as representative of any perturbation of the flux to the inorganic reservoir that carries the isotopic composition of organic matter.

Further simplifications follow by reexpressing the remaining terms in Eq. 11 using variables defined in the main text. First, recall that the mass m is decomposed into an unperturbed component m* and a perturbation Δm according toEmbedded Image(12)where, by definition,Embedded Image(13)Substitution of Eq. 12 into Eq. 9 and division of both sides by m* then yields, for Embedded ImageEmbedded Image(14)Recall also the definitions of the normalized mass perturbationEmbedded Image(15)the normalized input fluxEmbedded Image(16)and the turnover timeEmbedded Image(17)Inserting each of these expressions into Eq. 14, one obtainsEmbedded Image(18)The normalized mass perturbation evolves asEmbedded Image(19)which is obtained by taking the derivative of both sides of Eq. 13 and using Eqs. 15 to 17. Finally, the derivative of Eq. 4 in the main text provides the evolution of the normalized flux perturbationEmbedded Image(20)The differential equations 18 to 20 describe changes in the isotopic composition of marine inorganic carbon in response to the normalized flux J(t). The first term in parentheses in Eq. 18 encodes the tendency of δ1 to relax back to the initial value Embedded Image, and the second term describes the forcing of the system out of steady state by injection of carbon with isotopic composition Embedded Image. Because φ is unknown, the solution of this system requires not only three initial conditions but also one end condition:Embedded Image(21)Embedded Image(22)Embedded Image(23)Embedded Image(24)Equations 18 to 24 constitute the complete boundary value problem. The solution for φ, which is obtained numerically, provides the evolution of δ1(t) that matches its observed initial and final states while remaining consistent with Eqs. 18 to 20.

Turnover time, the organic burial fraction, and φc

Equation 17 shows that the turnover time τ0 equals the mass m* of inorganic carbon in the oceans divided by the input flux Embedded Image. In steady state, the latter is equal to the rate at which carbon is immobilized in marine sediments. The modern preindustrial reservoir of dissolved inorganic carbon is about 38,000 Pg (133). Sedimentary organic carbon and carbonate carbon are immobilized in sediments at the total preindustrial rate of approximately 0.28 Pg year−1 (134, 135). By dividing the reservoir size by the sequestration rate, one obtains τ0 ≃ 140 ky.

Krissansen-Totton et al. (21) have recently provided a comprehensive analysis of the evolution of the organic burial fraction through geologic time, commonly designated as f. Their study of the Phanerozoic portion of the record largely derives from an earlier analysis of Hayes et al. (43); both studies find f ≃ 0.23, which this paper assigns to φc. Hayes et al. (43) compute f in a series of contiguous intervals throughout the Phanerozoic, obtaining a sequence fi. The root mean square (rms) fluctuation of the fi’s is 0.037. The error assigned to φc is twice the rms fluctuation.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: This project was inspired by a great many discussions with S. Bowring. I thank him for generously sharing his insight into geologic time, the intellectual history of mass extinction, and the present understanding of the “sixth extinction.” I also thank T. Bosak, O. Devauchelle, R. Ferrari, L. Hatton, J. Hayes, P. Molnar, and R. Summons for helpful discussions. Funding: This work was supported by NASA Astrobiology grant NNA13AA90A and NSF grant EAR-1338810. Competing interests: The author declares that he has no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper. Additional data related to this paper may be requested from the author.