Asymmetry of extreme Cenozoic climate–carbon cycle events - Science Advances

jumianta.blogspot.com

INTRODUCTION

Paleoclimate proxy records reveal not only that Earth’s climate–carbon cycle system has changed substantially on time scales of many millions of years but also that it has experienced large, temporary disruptions on time scales of tens of thousands of years. Notable examples from within the Cenozoic [66 million years (Ma) to present] include the “hyperthermal” warming events of the early Eocene (1–13). This record of past climate–carbon cycle disruptions provides an observational window into the Earth system’s long-term response to anthropogenic forcing (14–16). However, important questions remain to be answered. To what extent does the record of past disruption reflect a proportionate response to external forcing, and to what extent does it reflect intrinsic self-amplification within the climate–carbon cycle system itself (17–19)? Because large disruptions appear to have been more common in some time periods than others (e.g., the Eocene), what are the general properties of the climate–carbon cycle system that determine the nature and magnitude of extreme events? The risk that a strongly nonlinear long-term Earth system response to anthropogenic forcing would pose for human civilization (20) adds urgency to these questions.

Cenozoic climate–carbon cycle fluctuations can be studied using δ¹⁸O and δ¹³C records from deep-sea benthic foraminifera (see Materials and Methods). Hyperthermal events are identified by paired negative excursions in δ¹⁸O and δ¹³C and have been interpreted as rapid global warming events caused by the release of isotopically depleted organic carbon into the surface environment. They further appear to have been paced by changes in the eccentricity of Earth’s orbit (4–11), although the precise mechanism is unclear. Proposed carbon sources for the largest hyperthermals include sedimentary methane hydrate (21, 22) or permafrost (23) reservoirs. However, many of the events likely reflect mechanisms of carbon release from relatively more accessible surficial reservoirs, such as dissolved organic carbon (24), that have persisted throughout much of the Cenozoic (8, 9).

The typical behavior of these hyperthermal disruptions is demonstrated in Fig. 1, which shows time series of benthic foraminiferal δ¹⁸O and δ¹³C from the early Eocene. The data are obtained from the global astronomically tuned Cenozoic composite record in (13). To isolate the sub-million-year fluctuations, we have subtracted a 1-Ma running mean. Finally, the empirical probability distribution of the fluctuations is also shown. Here, the hyperthermals manifest as extreme events in a probability distribution with an asymmetric non-Gaussian tail. The asymmetry quantifies a tendency toward negative excursions rather than positive excursions, suggesting that the climate–carbon cycle system exhibits a fundamental tendency toward extreme events involving global warming and oxidation of organic carbon. In this study, we quantify the evolution of this asymmetry throughout the Cenozoic and provide theoretical and analytical frameworks to explain the observed behavior.

Fig. 1 Climate–carbon cycle disruptions in the early Eocene, as recorded in benthic foraminiferal δ¹⁸O and δ¹³C data (13).

A 1-Ma running mean has been subtracted to isolate the sub-million-year fluctuations. (A and B) Time series and (C and D) histograms of the data points. The largest hyperthermals manifest as extreme events in an empirical probability distribution with an asymmetric non-Gaussian tail [near the asterisks in (C) and (D)]. This asymmetry quantifies an apparent bias toward extreme events involving global warming and oxidation of organic carbon. Note that the vertical axes decrease upward.

RESULTS

Cenozoic δ¹⁸O and δ¹³C fluctuations

Past studies of climate–carbon cycle disruptions have typically focused on individual, clearly identifiable “events.” Nevertheless, as Fig. 1 makes clear, there is a continuous spectrum of fluctuation sizes from these events all the way to the smallest fluctuations present in the data. Therefore, we use an alternative approach: studying the empirical probability distribution of all the available data points (as shown in Fig. 1, right). This is a widely used approach in the study of extreme weather and climate events on shorter time scales (25, 26) but has only rarely been applied to the study of paleoclimate proxy records (27, 28). In this context, it also has the additional advantage of being essentially insensitive to the specification of the underlying time scale.

We focus on robust features of empirical probability distributions that quantify the tendencies shown in Fig. 1. Letting X denote an arbitrary random variable, the asymmetry in the distribution p(X) can be characterized by the skewness$$S=E\left[{\left(\frac{X-E[X]}{\mathrm{\sigma }}\right)}^3\right]$$(1)where E denotes expectation and σ denotes the SD. The tendency toward extreme events can be characterized by the excess kurtosis (hereafter, kurtosis)$$K=E\left[{\left(\frac{X-E[X]}{\mathrm{\sigma }}\right)}^4\right]-3$$(2)

A positive kurtosis indicates that the probability distribution p(X) is heavy-tailed compared to the normal distribution.

Figure 2 shows the skewness and kurtosis of the δ¹⁸O and δ¹³C fluctuations in each epoch of the Cenozoic, together with 95% confidence intervals from a bootstrap analysis (see Materials and Methods). The Paleocene-Eocene Thermal Maximum has been removed because of its apparent uniqueness (9) and because its magnitude dwarfs the rest of the Eocene variability to the extent that it would hinder objective analysis of the more general behavior. Previous studies have considered a running skewness and kurtosis of portions of the Cenozoic δ¹⁸O record to quantify the nonsinusoidal nature of glaciation cycles (27, 28). Here, we choose to aggregate the data across epochs because we are focused on the large-scale trends. The skewness and kurtosis values for δ¹⁸O fluctuations in a given epoch should be very close to those of the temperature fluctuations in that epoch (with the sign of the skewness reversed; see Materials and Methods). The shorter-term variability in skewness and kurtosis throughout the Cenozoic would be interesting to explore in future work; nevertheless, it should not affect the results that we present here.

Fig. 2 Skewness and kurtosis of sub-million-year δ¹⁸O and δ¹³C fluctuations in the Cenozoic.

(A) Data organized by epoch (Pal = Paleocene, Eo = Eocene, Ol = Oligocene, Mio = Miocene, Plio = Pliocene, and Plei = Pleistocene). Error bars denote 95% confidence intervals (see Materials and Methods). All pre-Pliocene data exhibit negative skewness indicating an asymmetry that favors hyperthermal-like events. They also generally exhibit a positive kurtosis, indicating a greater tendency toward extreme events than would be expected from a normal distribution. (B) Data in skewness-kurtosis space. Shading indicates lower bounds for different classes of probability distributions: Distributions produced by correlated additive-multiplicative (CAM) noise processes cannot fall outside the white region (see Materials and Methods) (31), while unimodal distributions cannot fall in the dark gray region (30). Before the Pleistocene, the data tend to satisfy the more restrictive CAM bound. Many of the data are consistent with the lognormal distribution (black line), which is a further characteristic of multiplicative processes. These observations indicate that key dynamics of the system may be fruitfully described in terms of stochastic multiplicative noise.

Figure 2 reveals that δ¹⁸O and δ¹³C fluctuations have exhibited a substantial negative skewness and positive kurtosis throughout much of the Cenozoic. The negative skewness indicates an asymmetry favoring negative fluctuations of δ¹⁸O and δ¹³C, while the positive kurtosis indicates a greater tendency toward extreme events than would be expected from a normal distribution. These observations are not an expected consequence of orbital forcing (see Materials and Methods); we suggest instead that they arise from intrinsic features of the climate–carbon cycle system. They quantify the bias toward hyperthermal-like extreme events observed in Fig. 1; the fact that this bias is not unique to the Eocene is in line with previous suggestions that Eocene hyperthermal events reflected mechanisms persisting throughout much of the Cenozoic (8, 9).

Although the skewness of the δ¹⁸O and δ¹³C fluctuations varies in magnitude over time, its negative sign persists throughout all epochs before the Pliocene (5.3 Ma to 2.6 Ma). During the Pliocene, the δ¹⁸O fluctuations instead become positively skewed. This change in sign is suggestive of a “switch” in the coupling of the climate and the carbon cycle, perhaps related to the onset of Northern Hemisphere glaciation (29). In the Pleistocene (2.6 Ma to present), the kurtosis of both δ¹⁸O and δ¹³C fluctuations decreases substantially; this indicates a lessened susceptibility to extreme events and may thus reflect an increase in the stability of the climate–carbon cycle system.

These observations become more intriguing when one considers the predicted skewness-kurtosis relationships for different classes of probability distributions. For example, Klaassen et al. (30) have shown that all unimodal distributions must satisfy$$K\ge S^2-\frac{186}{125}$$(3)

A much more restrictive bound exists for the distribution of fluctuations produced by stochastic processes involving correlated additive-multiplicative (CAM) noise (discussed further below) (26, 31, 32)$$K\ge \frac{3}{2}{S}^2-r$$(4)where r = 0 for single-variable systems and r has a small positive value for systems with multiple variables (see Materials and Methods) (32).

Figure 2 shows that the δ¹⁸O and δ¹³C fluctuations before the Pleistocene not only satisfy the unimodal bound (Eq. 3) but also tend to satisfy the much more restrictive one-variable CAM bound (Eq. 4). Furthermore, many of the data points are consistent with the lognormal distribution (see Materials and Methods), which emerges generally from a range of multiplicative processes (33). These observations suggest that key dynamics of the climate–carbon cycle system may be fruitfully described in terms of stochastic multiplicative noise.

Multiplicative noise in the climate–carbon cycle system

Stochastic models were first applied to the study of climate variability by Hasselmann (34), who used a model of the form$$\frac{\mathit{dx}}{\mathit{dt}}=-\frac{1}{\mathrm{\tau }}x+\mathrm{\nu }\mathrm{\eta }(t)$$(5)to understand the “red” power spectrum of many weather and climate time series. Here, x represents the variable of interest, τ is the time scale on which negative feedbacks tend to return the system toward x = 0, and η(t) is Gaussian white noise (see Materials and Methods). Note that the white noise term does not represent “true” stochasticity but is rather an approximation of the combined effects of many deterministic fluctuations that decorrelate on a time scale much shorter than that of the long-term climate variations being considered.

In Eq. 5, the intrinsic fluctuations νη(t) have an amplitude that is independent of the system state x. This is referred to as “additive noise” and causes the probability distribution of the output fluctuations, p(x), to be Gaussian (34, 35). In contrast, if the amplitude of the intrinsic fluctuations depends on the system state, we obtain the simple “multiplicative noise” model$$\frac{\mathit{dx}}{\mathit{dt}}=-\frac{1}{\mathrm{\tau }}x+f(x)\mathrm{\eta }(t)$$(6)

If f(x) is an increasing function of x, the influence of the noise is greater when x is larger. Consequently, Eq. 6 generates probability distributions p(x) that are asymmetric and have heavier tails than Gaussian distributions. Models that include multiplicative noise have been previously applied to study a wide range of climate problems (36–39) but on much faster time scales than those we consider here.

A useful special case of a multiplicative noise model is obtained by linearizing the state dependence in Eq. 6 around x = 0. This yields a simple one-variable CAM noise model$$\frac{\mathit{dx}}{\mathit{dt}}=-\frac{1}{\mathrm{\tau }}x+\mathrm{\nu }(x+c)\mathrm{\eta }(t)$$(7)

Linear CAM noise models have been used to study extreme weather events (26, 31, 32) and are attractive, in part, due to their analytical tractability: The steady-state probability distribution p(x), as well as the kurtosis-skewness lower bound (4), can be straightforwardly derived (see Materials and Methods). While the steady-state probability distribution for Eq. 5 is Gaussian, the steady-state distribution for Eq. 7 has an asymmetric non-Gaussian tail, in agreement with Fig. 1. It further has kurtosis and skewness values that satisfy $K\ge \frac{3}{2}{S}^2$, consistent with the general behavior of δ¹⁸O and δ¹³C fluctuations before the Pleistocene (Fig. 2).

While the CAM noise model predicts a lower bound (Eq. 4), an exact kurtosis-skewness relationship emerges in the case of pure multiplicative fluctuations on a time scale much faster than the damping time scale τ. In this case, Eq. 7 reduces to$$\frac{\mathit{dx}}{\mathit{dt}}=\mathrm{\nu }x\mathrm{\eta }(t)$$(8)and x will be lognormally distributed (see Materials and Methods). The fact that many of the data points in Fig. 2 are consistent with lognormal behavior thus underscores the potential significance of multiplicative noise in generating the hyperthermal-like extreme events observed throughout the Cenozoic. A schematic summarizing the behaviors of these simple models is shown in Fig. 3.

Fig. 3 Schematic summarizing stochastic models discussed in the text.

(A) Additive noise model (34): Here, the noise amplitude is independent of the system state. This produces a Gaussian distribution of fluctuations, with K, S = 0. (B) CAM noise (Eq. 7) generates asymmetric non-Gaussian distributions satisfying $K\ge \frac{3}{2}{S}^2$, consistent with the pre-Pleistocene data in Fig. 2. (C) Undamped multiplicative noise produces a lognormal distribution of fluctuations, which also has an asymmetric non-Gaussian tail. There exists an exact parametrizable K(S) relationship (see Materials and Methods): It is plotted in Fig. 2 and intersects a number of the data points.

Multiplicative noise could replicate the pre-Pliocene asymmetry favoring hyperthermal-like events if the amplitude of the intrinsic fluctuations increases as the δ¹⁸O and δ¹³C anomalies become more negative. What could be responsible for such a relationship in the global climate–carbon cycle system? One attractive possibility is that it reflects the effects of temperature (which is inversely related to δ¹⁸O) on biological and chemical reaction rates (see Materials and Methods). The fast deterministic fluctuations that are approximated as intrinsic random noise involve biological and chemical processes, whose rates increase with temperature. In the context of Eq. 6, increasing these rates corresponds to increasing the amplitude of the random noise term; thus, it seems reasonable that increased global temperatures could increase the amplitude of intrinsic fluctuations in the climate–carbon cycle system.

If the amplitude of intrinsic fluctuations in the climate–carbon cycle system exhibits a positive correlation with global temperature before the Pliocene, this should have left additional signatures in the geochemical record. We investigate this by dividing the Cenozoic δ¹⁸O time series into 0.5-Ma bins and testing for a negative relationship between the mean δ¹⁸O and an estimate of the amplitude of the intrinsic δ¹⁸O fluctuations in each bin (see Materials and Methods). For each epoch, we compute rank correlation coefficients between these two variables, together with significance levels from a Monte Carlo permutation test. Table 1 shows that a negative relationship between the mean δ¹⁸O and the amplitude of the intrinsic fluctuations is exhibited in each epoch. For the Eocene and Miocene, this relationship is statistically significant with P < 0.05, and combined P values across all four epochs are also significant (see Materials and Methods).

Table 1 Relationship between mean δ¹⁸O and amplitude of intrinsic fluctuations.

This is quantified in terms of a rank correlation (see Materials and Methods). In each epoch, the amplitude of underlying fluctuations increases with decreasing δ¹⁸O, consistent with the multiplicative noise hypothesis. For the Eocene and Miocene, this relationship is statistically significant when considering only the data from that epoch (P < 0.05), but combined P values across all four epochs are also statistically significant (see Materials and Methods).

View this table:

These results suggest that the amplitude of intrinsic fluctuations in the climate–carbon cycle system may have increased with temperature before the Pliocene, consistent with the multiplicative noise hypothesis stated above. Even though changes in global ice volume make an important contribution to the δ¹⁸O signal starting in the Oligocene, this inference remains robust due to our use of rank correlations. Although the presence or absence of ice sheets modifies the δ¹⁸O − T relationship, to a good approximation, the relationship remains linear within each epoch (see Materials and Methods). Then, regardless of the details, a negative rank correlation in Table 1 will correspond to a positive rank correlation between mean temperature and the amplitude of fluctuations with precisely the same magnitude and significance level (see Materials and Methods).

Stochastic climate–carbon cycle model

To better understand how temperature-driven multiplicative noise in the climate–carbon cycle system would generate asymmetric extreme events in δ¹⁸O and δ¹³C like those in Fig. 1, we develop a simple stochastic numerical climate–carbon cycle model (see Materials and Methods). The model considers the evolution of global surface temperature and surficial inorganic carbon, aspects of ocean chemistry, the CO₂ greenhouse effect, and the long-term weathering feedback, producing δ¹⁸O and δ¹³C output time series. It assumes that the amplitude of intrinsic random fluctuations in the surficial inorganic carbon reservoir, which are driven by temporary imbalances in the global production and oxidation of organic carbon, increases as the global mean temperature increases. We further include stochastic fluctuations in global mean temperature and assume that they are partially correlated with these carbon cycle fluctuations.

Figure 4 shows the result of forcing this model with 400-thousand-year (ka) eccentricity variations (see Materials and Methods) using an ensemble of 100 trajectories; a single trajectory has been highlighted in black. The model generates a full spectrum of variability with hyperthermal-like extreme events (paired negative δ¹⁸O and δ¹³C excursions) that have a tendency to occur near eccentricity maxima, consistent with observations (4–11). Each individual extreme event occurs due to release of isotopically depleted organic carbon into the ocean-atmosphere system, consistent with previous suggestions (5–13, 21–23); the tendency for orbital pacing arises through the effect of the eccentricity forcing on the noise amplitude (see Materials and Methods). Across the ensemble of different model realizations, key behaviors of the pre-Pliocene data in Fig. 2 are reproduced: The distribution of fluctuations in both proxies is negatively skewed, kurtosis tends to be positive, the skewness and kurtosis values behave in accordance with the lower bound for multivariable CAM noise models (Eq. 4; see Materials and Methods), and many of the data fall near the lognormal line. Finally, the average slope of δ¹⁸O fluctuations with respect to the δ¹³C fluctuations is also consistent with the pre-Pliocene data (see Materials and Methods).

Fig. 4 Numerical model results.

(A) Ensemble of 100 trajectories obtained by forcing the stochastic climate–carbon cycle model with 400-ka eccentricity variations. A single trajectory is highlighted in black. The model generates a full spectrum of variability with hyperthermal-like extreme events (paired negative δ¹⁸O and δ¹³C excursions) that have a tendency to occur near eccentricity maxima, consistent with observations (4–11). (B) Skewness and kurtosis values for the different ensemble trajectories, together with the bound for unimodal probability density functions (30), the multivariable CAM noise bound (Eq. 4; see Materials and Methods), and the relationship for the lognormal distribution; we observe similar behavior to the pre-Pliocene observations in Fig. 2. (C) Scatterplot of δ¹⁸O versus δ¹³C from all of the model output together with a linear regression and the corresponding regression lines for the pre-Pliocene observations (see Materials and Methods); there is again good agreement.

DISCUSSION

In this work, we have quantified general trends in the behavior of extreme events in the climate–carbon cycle system throughout the Cenozoic. We found that sub-million-year fluctuations in epochs before the Pliocene exhibited a fundamental asymmetry, favoring extreme events involving negative excursions in δ¹⁸O and δ¹³C. This is consistent with an implicit existing understanding that extreme climate–carbon cycle events have generally been hyperthermal-like. The fluctuations also tended to exhibit positive kurtosis, indicating an amplification of extreme events (i.e., a heavier tail) relative to the normal distribution. The quantitative persistence of both behaviors throughout much of the Cenozoic, as shown in Fig. 2, suggests that hyperthermal-like disruptions arise not only as interesting individual events but also as a general consequence of intrinsic features of the climate–carbon cycle system.

Our results show that the behavior of extreme climate–carbon cycle events throughout the Cenozoic is well described in terms of stochastic multiplicative noise. Stochastic multiplicative noise fundamentally generates asymmetric non-Gaussian fluctuations, in quantitative agreement with the observations. For example, the skewness and kurtosis of the observed Cenozoic δ¹⁸O and δ¹³C fluctuations tend to satisfy the lower bound $K\ge \frac{3}{2}{S}^2$ for fluctuations produced by CAM noise (Fig. 2), which is much more restrictive than the requirement for unimodal distributions. Furthermore, intrinsic climate–carbon cycle fluctuations appear to increase in amplitude with decreasing δ¹⁸O before the Pliocene, exactly as expected for multiplicative noise (Table 1). Last, a numerical climate–carbon cycle model in which the amplitude of fluctuations in the surficial carbon inventory increases with temperature is able to reproduce asymmetric hyperthermal-like extreme events, observed skewness-kurtosis relationships, δ¹⁸O − δ¹³C slopes, and the observed pacing of hyperthermal-like events by changes in orbital parameters (Fig. 4).

Beyond reproducing observations, the multiplicative noise perspective likely offers fundamental insight into the real climate–carbon cycle system. Past modeling work has focused on understanding how carbon may be released from buried sedimentary sources (21–23) and on deducing the nature of the carbon release events responsible for specific isotopic excursions (15, 40). Multiplicative noise, on the other hand, provides a dynamical explanation of how and why hyperthermal-like events throughout the Cenozoic could have arisen generally from processes of carbon redistribution between Earth’s relatively accessible surficial reservoirs. Specifically, it suggests that fluctuating imbalances in the global production and oxidation of organic carbon were amplified in the direction of carbon release by multiplicative effects, potentially due to the temperature dependence of biological and chemical reaction rates. The lognormal-like behavior of many observed (Fig. 2) and simulated (Fig. 4) data then further indicates that those multiplicative bursts were largely underdamped with respect to long-term stabilizing weathering feedbacks, consistent with a substantial time scale separation of the underlying processes.

This study also provides a new framework within which to investigate differences between the different epochs of the Cenozoic. What is the origin of the many different behaviors observed in Fig. 2? For example, why do δ¹³C fluctuations in the Eocene and Miocene exhibit a more negative skewness and a greater kurtosis than the corresponding δ¹⁸O fluctuations while this trend is reversed in the Paleocene, and why is the magnitude of both the skewness and kurtosis lower in the Oligocene? The Pliocene fluctuations appear consistent with multiplicative noise, but the changed sign of the δ¹⁸O asymmetry remains to be addressed: Is this a consequence of the onset of Northern Hemisphere glaciation (29)? On the other hand, the much lower kurtosis of the Pleistocene system is inconsistent with multiplicative noise, suggesting that it has been in some way more stable. The development of glacial cycle oscillations (41–43) may have “seized control” of the climate–carbon cycle system, damping the processes that earlier led to the asymmetric amplification of extreme events. This suggests that this asymmetric amplification may return as anthropogenic warming continues and the Northern Hemisphere ice sheets disappear, making the Earth system more susceptible to extreme warming events occurring on time scales of tens of thousands of years.

MATERIALS AND METHODS

δ¹⁸O and δ¹³C fluctuations by epoch

In this study, we use δ¹⁸O and δ¹³C data from the Cenozoic Global Reference benthic foraminifer carbon and oxygen Isotope Dataset (CENOGRID) (13). δ¹⁸O is inversely related to deep-sea temperature (see below), and δ¹³C records changes in the carbon cycle. We isolate sub-million-year fluctuations by subtracting a 1-Ma moving average from the data. For each geologic epoch within the Cenozoic (44), we calculate the skewness and kurtosis of the empirical distribution of fluctuations: Because the number of data points per epoch is reasonably large (>2000), we interpret the expected values in Eqs. 1 and 2 as straightforward sample averages. For non-Gaussian fluctuations, the sampled skewness and kurtosis across a given interval is strongly affected by how many of the more extreme events occur within that interval. Dividing the data by epoch allows us to keep these intervals as large as possible while still capturing long-term trends that can be related to the important changes occurring between epoch boundaries (e.g., the onset of Northern Hemisphere glaciation in the Pliocene).

We obtain 95% confidence intervals for our skewness and kurtosis estimates using a bootstrap method: Letting N denote the number of data points in a given epoch, we create bootstrap samples of size N by randomly sampling from the observations with replacement, and then calculate that sample’s skewness and kurtosis. Repeating this procedure 1000 times yields approximate error distributions for skewness and kurtosis values: Denoting the statistic of interest as x, the 95% confidence interval is [c₁, c₂], where P(x < c₁) = P(x > c₂) = 0.025.

δ¹⁸O, temperature, and ice volume

The relationship between temperature and the isotopic composition of foraminiferal calcite is typically parametrized as (45)$$T=a-b({\mathrm{\delta }}^{18}{\mathrm{O}}_{\text{calcite}}-{\mathrm{\delta }}^{18}{\mathrm{O}}_{\text{water}})$$(9)where δ¹⁸O_calcite and δ¹⁸O_water are the isotopic compositions of the calcite and the surrounding water, respectively, and a and b are constants. Because the growth of ice sheets increases δ¹⁸O_water, the benthic δ¹⁸O signal reflects both changes in temperature and in global ice volume. The relative importance of each factor changes throughout the Cenozoic notably with the onset of Southern Hemisphere glaciation at the start of the Oligocene and the onset of Northern Hemisphere glaciation at the start of the Pliocene. Nevertheless, previous work suggests that the expression$$T=\mathrm{\alpha }-\mathrm{\beta }{\mathrm{\delta }}^{18}{\mathrm{O}}_{\text{calcite}}$$(10)remains a good approximation: Within different epochs, the presence or absence of ice sheets modifies the slope β and the offset α (46).

Our analysis of the skewness and kurtosis of δ¹⁸O fluctuations in Fig. 2 stands independently of the δ¹⁸O − T relationship. The linear relationship (Eq. 10), however, greatly aids physical interpretation. As long as α and β can be approximated as constant within a given epoch, the fluctuations in T have a skewness and kurtosis of precisely the same magnitude as the fluctuations of δ¹⁸O (the skewness will have the opposite sign).

Role of orbital forcing

On long time scales, the climate–carbon cycle system is forced by quasiperiodic variations in Earth’s orbital parameters. These variations have been calculated in detail (47, 48), and their imprint is evident in the δ¹⁸O and δ¹³C records (13). Precisely how the orbital variations actually force the climate–carbon cycle system has not yet been settled; previous studies have highlighted the likely importance of low- to mid-latitude insolation changes (4, 49).

We evaluate whether the orbital forcing could be responsible for the asymmetry and the non-Gaussian tails in δ¹⁸O and δ¹³C fluctuations (Fig. 2) by analyzing the statistics of the orbital solutions calculated in (47). We would consider the orbital forcing to be “responsible” for these observations if the observations can be explained by a simple linear response of the climate–carbon cycle system without nonlinear amplification. We generate the time series of insolation from 100 Ma to present sampled at a 1-ka time step; because we are focused on statistics, our argument does not depend on the specific time intervals chosen. The insolation variations at a given latitude are very close to Gaussian: This is demonstrated in fig. S1 for insolation at the equator and at 45^∘N. There is a very modest skewness (much smaller than the skewnesses in Fig. 2), and a negative kurtosis, rendering these variations insufficient for explaining the observation of a substantial skewness and positive kurtosis in the δ¹⁸O and δ¹³C fluctuations in terms of a linear response.

Note that the average insolation received by Earth over a whole year does exhibit fluctuations with a substantial skewness and kurtosis; a histogram is plotted in fig. S2. This is a straightforward consequence of near-Gaussian fluctuations in the eccentricity, e, and the mean annual insolation scaling as $1/\sqrt{1-e^2}$ (50). However, this behavior is also insufficient for explaining the observed asymmetry and heavy tails in the fluctuations in δ¹⁸O and δ¹³C for multiple reasons. First, when considering orbital forcing of the climate–carbon cycle system, the mean annual insolation is likely not the relevant quantity, as discussed above. Second, the magnitude of these variations is far too small (∼0.8 W/m²) to account for the magnitude of the observed extreme events (e.g., in δ¹⁸O) without nonlinear amplification of some kind. Third, the skewed heavy tail in the mean insolation represents eccentricity variations on time scales ≳100 ka, while the skewed non-Gaussian tail in the δ¹⁸O and δ¹³C observations represents fluctuation events occurring on shorter time scales. Last, the kurtosis of the mean annual insolation variations falls far below that predicted by the CAM bound $K\ge \frac{3}{2}{S}^2$ (fig. S2): Even without the problems discussed above, it would still need to be explained why the observations behave differently. These considerations suggest that mechanisms intrinsic to the climate–carbon cycle system play a dominant role in generating the observed asymmetric non-Gaussian tails in the δ¹⁸O and δ¹³C fluctuations.

Stochastic multiplicative noise theory

In Eqs. 5 to 8, η(t) is delta-correlated Gaussian white noise satisfying 〈η(t₁)η(t₂)〉 = δ(t₁ − t₂) and 〈η(t)〉 = 0. Note also that, throughout this paper, we have chosen to interpret stochastic differential equations using the Itô calculus; conversion to the related Stratonovich calculus is straightforward (35).

The steady-state probability distribution for the additive noise model (Eq. 5) is straightforwardly obtained by integrating the corresponding Fokker-Planck equation (35): It is the Gaussian distribution$$p(x)=\frac{1}{\sqrt{2\mathrm{\pi }{\mathrm{\sigma }}^2}}\text{exp} (-\frac{{(x-\mathrm{\mu })}^2}{2{\mathrm{\sigma }}^2})$$(11)with μ = 0 and σ² = τν²/2. The steady-state distribution for the one-variable CAM model (Eq. 7) is obtained similarly, yielding$$p(x)\propto \text{exp} (-\frac{2c}{\mathrm{\tau }{\mathrm{\nu }}^2(x+c)}){(x+c)}^{-2(1+\frac{1}{\mathrm{\tau }{\mathrm{\nu }}^2})}$$(12)

The $K\ge \frac{3}{2}{S}^2$ relationship (Eq. 4), as well as the steady-state distribution (Eq. 12), is derived in (32). Because of the importance of these results to this paper and because we have used slightly different notation as well as a different stochastic calculus, the Supplementary Materials include a derivation of both results directly from Eq. 7. The relationship $K\ge \frac{3}{2}{S}^2-r$ can be further obtained for multivariable linear systems with CAM noise under the assumption that the operator describing the deterministic evolution is nonnormal; for a derivation, and a discussion of the validity of this assumption in geophysical contexts, the reader is referred to (32).

The lognormal distribution arises from a range of multiplicative processes, in part due to the central limit theorem (33). In the context of Eq. 8, its appearance can be understood by substituting y = log x and noting that y evolves according to an additive noise process, thus obeying the normal distribution (35, 51). The solution to Eq. 8 obeys the lognormal distribution$$p(x)=\frac{1}{\mathrm{\sigma }x\sqrt{2\mathrm{\pi }}}\text{exp} (-\frac{{(\text{ln} x-\mathrm{\mu })}^2}{2{\mathrm{\sigma }}^2})$$(13)where $\mathrm{\mu }=\text{log} (x(t=0))-\frac{1}{2}{\mathrm{\nu }}^{2}t$ and σ² = ν²t. The kurtosis-skewness relationship can then be described parametrically through the expressions (33)$$S=(\text{exp} ({\mathrm{\sigma }}^2)+2)\sqrt{\text{exp} ({\mathrm{\sigma }}^2)-1}$$(14)$$K=\text{exp} (4{\mathrm{\sigma }}^2)+2\text{exp} (3{\mathrm{\sigma }}^2)+3\text{exp} (2{\mathrm{\sigma }}^2)-6$$(15)

The plotted line in Fig. 2 incorporates both the cases where log X and log (− X) are normally distributed. In the latter case, the skewness (Eq. 14) changes sign.

Effect of temperature on reaction rates

Multiplicative noise could replicate the pre-Pliocene asymmetry favoring hyperthermal-like events if the amplitude of intrinsic fluctuations increases as the δ¹⁸O anomaly decreases (i.e., temperature increases). The deterministic processes in the climate–carbon cycle system that we are approximating as random noise involve biological and chemical processes, whose rates would increase as temperature increases. The rates of many chemical reactions increase with temperature according to the Arrhenius relationship (52)$$k\propto \text{exp} (-\frac{E_{\mathrm{a}}}{k_{\mathrm{b}}T})$$(16)where E_a is an activation energy and k_b is Boltzmann’s constant. Similar behavior may apply to the biologically mediated reactions that constitute the global carbon cycle (53, 54). We therefore argue that it is reasonable to expect the amplitude of intrinsic fluctuations within the global carbon cycle to increase with temperature; our analysis of the δ¹⁸O record provides further tentative evidence supporting this (Table 1).

Further signatures of multiplicative noise in Cenozoic δ¹⁸O data

To interrogate the observations for further signatures of multiplicative noise, we investigate how the amplitude of the intrinsic fluctuations in the δ¹⁸O record changes with the 1-Ma mean of δ¹⁸O in all of the data before the Pliocene. The simplest possible metric of this amplitude would be the SD of δ¹⁸O about the long-term mean, but, across any given time interval, this will be strongly affected by the number of extreme events that occur. Because these extreme events almost uniformly occur in the direction of negative δ¹⁸O, we can remove them from our estimate of the magnitude of the intrinsic fluctuations by considering only the positive fluctuations above the mean; a similar approach was used in (49). We divide the δ¹⁸O time series into 0.5-Ma bins and, for each bin, calculate the mean δ¹⁸O and the SD of the positive fluctuations.

We test for a monotonic relationship between the binned means and fluctuation amplitudes across each epoch by calculating Spearman rank correlation coefficients (55). Significance levels are calculated using a Monte Carlo permutation test: randomly reorder the relationships between the binned means and amplitudes, and then recalculate the correlation coefficients. Repeating this procedure 10,000 times yields a distribution of rank correlation coefficients under the null hypothesis that the mean δ¹⁸O and the fluctuation amplitude in each bin are uncorrelated. The significance levels for the observed rank correlations are then straightforwardly calculated from this distribution.

We find a negative relationship between δ¹⁸O and the fluctuation amplitude, consistent with the behavior required to generate the asymmetric non-Gaussian tails in Fig. 2. For the Eocene and Miocene epochs, these negative rank correlations are statistically significant with P < 0.05. Although the negative δ¹⁸O fluctuation relationships observed in the Paleocene and Oligocene epochs are not statistically significant at this level given those data alone, we note that combined P values that take into account all four of the epochs considered are very small: P < 9 × 10⁻⁶ using Fisher’s method (56) and P < 4 × 10⁻⁵ using the harmonic mean (57).

Because of the negative relationship between δ¹⁸O and temperature, this result is consistent with the temperature-driven multiplicative noise hypothesis. While decreasing δ¹⁸O corresponds to increasing T throughout the Cenozoic, the precise shape of this relationship has been affected by the presence of ice sheets, starting in the Oligocene. Nevertheless, the use of a rank correlation means that the results in Table 1 can be robustly interpreted in terms of temperature. As long as δ¹⁸O(T) is monotonically decreasing within each epoch considered, the rank orders of the binned mean values will stay the same. As long as δ¹⁸O(T) is approximately linear within each epoch considered [as suggested, e.g., in (46)], the rank orders of the binned fluctuation amplitudes will stay the same. If the rank orders remain the same, the negative correlation coefficients for δ¹⁸O in Table 1 become positive correlation coefficients for T with precisely the same magnitude and significance levels.

Last, it is possible that the Miocene result in Table 1 is affected by the near–ice-free conditions of the mid-Miocene climatic optimum: The slope of the δ¹⁸O − T relationship could have changed during this time. Because the δ¹⁸O − T relationship becomes less steep in an ice-free period, this could introduce a positive bias into the correlation between mean δ¹⁸O and δ¹⁸O fluctuations. However, since we have observed a negative correlation (Table 1), our basic result (amplitude of fluctuations increasing with global temperature) remains robust.

Stochastic climate–carbon cycle model

Our stochastic climate–carbon cycle model considers the total amount of ocean-atmosphere inorganic carbon, I, the deviation of the global mean surface temperature from a long-term stable state, ΔT, and the amount of ocean-atmosphere inorganic ¹³C, I₁₃. Note that “long-term” here refers to time scales of millions of years or greater. We do not consider changes in this long-term stable state (e.g., due to tectonic processes), as we are focused on the sub-million-year fluctuations.

On time scales of hundreds of thousands of years, I is widely thought to be controlled by a stabilizing feedback provided by the weathering of carbonate and silicate rocks (58, 59). Defining I₀ as the long-term steady-state value of I (all parameter values are given in the next section), this stabilizing feedback can be simply parametrized as$$\frac{\mathit{dI}}{\mathit{dt}}=-\frac{(I-I_0)}{\mathrm{\tau }}$$(17)where τ is the characteristic time scale of the weathering feedback. Following the analysis in the main text, we include fluctuations in I that arise from imbalances in the production and oxidation of organic carbon and assume that the amplitude of these fluctuations increases with temperature. We parametrize this as CAM noise, leading to the equation$$\frac{\mathit{dI}}{d\mathrm{t}}=-\frac{(I-I_0)}{\mathrm{\tau }}+{\mathrm{\nu }}_C(\mathrm{\Delta }T+c)\mathrm{\eta }(t)$$(18)

Here, η(t) is a Gaussian white noise process as described above, and ν_C and c control the strength of the temperature dependence and the amplitude of the noise at ΔT = 0.

The global reservoir of organic carbon, which grows when η(t) < 0 (net production) and shrinks when η(t) > 0 (net oxidation), is left implicit. We consider this global reservoir to consist of the sum total of relatively accessible surficial organic carbon stocks, such as dissolved organic carbon (8, 24). This implicit formulation is reasonable, in part, because the fluctuation term ν_C(ΔT + c)η(t) in Eq. 18 has mean zero and does not contribute to any mean drift in I; in other words, on average, it acts as neither a source nor sink of inorganic carbon.

The deterministic evolution of global mean surface temperature is determined by the balance of incoming and outgoing radiation; because outgoing radiation is approximately linear in surface temperature and log CO₂ for a wide range of parameters (60), this can be parametrized as$$\frac{d\mathrm{\Delta }T}{\mathit{dt}}=\frac{1}{C}(-a_1(\mathrm{\Delta }T)+a_2\text{log} (\frac{P(I)}{P_0}\left)\right)$$(19)

Here, C denotes the surface heat capacity, P denotes the atmospheric CO₂ concentration, P₀ is the steady-state CO₂ concentration, and a₁ and a₂ are constants. P is obtained directly from I and ocean carbonate chemistry under the assumption that total alkalinity remains constant (see the Supplementary Materials). We also introduce stochastic fluctuations in global mean surface temperature that are partially correlated with those in Eq. 18. The full temperature evolution equation is then$$\frac{d\mathrm{\Delta }T}{\mathit{dt}}=\frac{1}{C}(-a_1(\mathrm{\Delta }T)+a_2\text{log} (\frac{P(I)}{P_0}\left)\right)+{\mathrm{\nu }}_T(\mathrm{\Delta }T+c)\mathrm{\eta }(t)+\mathrm{\mu }\mathrm{\xi }(t)$$(20)

Here, η(t) is the same Gaussian white noise process as in Eq. 18, while ξ(t) is Gaussian white noise independent of η(t). Their amplitudes are controlled by the parameters ν_T and μ.

Our model for the evolution of δ¹³C values follows in spirit from those in (24, 61), although here we express our equations in terms of an explicit ¹³C variable for reasons of numerical stability. The evolution of ocean-atmosphere inorganic ¹³C, I₁₃, follows mechanistically from Eq. 18. We decompose the weathering feedback term −(I − I₀)/τ into an incoming flux I₀/τ of carbonate carbon with isotopic composition δ_c and an outgoing flux I/τ with the isotopic composition of the surficial inorganic carbon reservoir. The deterministic evolution of I₁₃ is then given by$$\frac{d{I}_{13}}{\mathit{dt}}=-\frac{(I_{13}-I_0{R}_{\mathrm{c}})}{\mathrm{\tau }}$$(21)where R_c represents the ¹³C/(¹²C+¹³C) ratio corresponding to δ_c. Neglecting the small difference between ¹³C/¹²C and ¹³C/(¹²C+¹³C), this conversion is carried out using$$R=R_{\text{std}}(1+\frac{\mathrm{\delta }}{1000})$$(22)where R_std represents the Vienna Pee Dee Belemnite (VPDB) standard.

Last, Eq. 21 also needs to account for the stochastic fluctuations in Eq. 18, ν_C(ΔT + c)η(t). Letting δ_i denote the isotopic composition of the inorganic carbon reservoir, these fluctuations would either remove carbon with an isotopic composition δ_i − ε (where ε > 0 denotes fractionation) or add organic carbon with an isotopic composition δ_o. On the sub-million-year time scales and for the relatively small changes we are concerned with, it is reasonable to assume that δ_o = δ_i − ε, where ε is a constant, leading directly to$$\frac{d{I}_{13}}{\mathit{dt}}=-\frac{(I_{13}-I_0{R}_{\mathrm{c}})}{\mathrm{\tau }}+{\mathrm{\nu }}_{\mathrm{C}}(\mathrm{\Delta }T+c)R_{\mathrm{o}}\mathrm{\eta }(t)$$(23)where R_o denotes the ¹³C/(¹²C+¹³C) ratio corresponding to δ_o.

The model is fully specified by Eqs. 18, 20, and 23. Once it has been run, a δ¹⁸O time series is obtained (with an arbitrary offset) as$${\mathrm{\delta }}^{18}\mathrm{O}(t)=-\frac{\mathrm{\Delta }T(t)}{4.8}$$(24)where the conversion constant (45) is for ice-free conditions (e.g., the Eocene). The appropriate linear conversion constants for different time periods within the Cenozoic can be found in (46). As noted earlier, as long as the relationship remains linear, the choice of conversion constant does not affect the skewness and kurtosis of the empirical distribution of fluctuations. A δ¹³C time series is obtained as$${\mathrm{\delta }}^{13}\mathrm{C}(t)=(\frac{I_{13}(t)/I(\mathrm{t})}{R_{\text{std}}}-1)\times 1000$$(25)

The model is implemented in Julia using the package DifferentialEquations.jl (62) and integrated using an Euler-Maruyama algorithm. Full model code is available at

https://github.com/arnscheidt/asymmetric-cenozoic-extreme-events.

Model parameter values

The parameter values used in the model are I₀ = 38,000 Pg, τ = 100 ka, C = 2 × 10⁻⁸ J m⁻² K⁻¹, P₀ = 400 μatm, δ_c = 1‰, δ_o = −25‰, ν_T = 0.2 year^−1/2, ν_C = 1.0 year^−1/2, c = 1.0 K, μ = 0.4 K/year^−1/2, a₁ = 2.2 W m⁻² K⁻¹. Because the steady state of the deterministic temperature evolution equation (Eq. 19) reduces to$$\mathrm{\Delta }T=\frac{a_2}{a_1}\text{log} \left(\frac{P}{P_0}\right)$$(26)it is convenient to let a₂ be expressed in terms of a₁ and the long-term temperature response of the Earth system to a doubling of CO₂, λ$$\frac{a_2}{a_1}=\frac{\mathrm{\lambda }}{\text{log}(2)}$$(27)

Here, we use λ = 5 K, consistent with its interpretation as an “Earth system sensitivity,” e.g., in the sense of (63).

Forcing the model with changes in insolation

Periodic insolation forcing is implemented by modifying the temperature evolution equation to read$$\frac{d\mathrm{\Delta }T}{\mathit{dt}}=\frac{1}{C}(-a_1(\mathrm{\Delta }T-F(t))+a_2\text{log}(\frac{P(I)}{P_0}\left)\right)+{\mathrm{\nu }}_{\mathrm{T}}(\mathrm{\Delta }T+c)\mathrm{\eta }(t)+\mathrm{\mu }\mathrm{\xi }(t)$$(28)where F(t) is a time-varying function that sets the “effective steady-state” temperature. With this formulation, for I near I₀, the system will radiatively adjust toward the temperature of T₀ + F(t). For the demonstration in Fig. 4, we use a simple eccentricity-like forcing$$F(t)=3\text{sin} \left(\frac{2\mathrm{\pi }t}{400 \text{ka}}\right)$$(29)

The tendency of the extreme events in Fig. 4 to be paced by variations in eccentricity arises because the amplitude of the intrinsic random fluctuations in our climate–carbon cycle model increases at higher temperatures, making extreme events more likely. Note that the annual mean insolation variations due to eccentricity are not directly large enough to cause surface temperature changes of the magnitude implied by Eq. 29. However, the climate–carbon cycle system likely contains mechanisms that transfer power from precession to eccentricity frequencies (49, 64). If hyperthermal-like extreme events occur due to multiplicative noise in the climate–carbon cycle system, then as long as eccentricity interacts in some way with the noise amplitude, the extreme events will tend to be paced by it.

Kurtosis and skewness of model output trajectories

The skewness and kurtosis values plotted in Fig. 4 obey the bound in Eq. 4 with r ≃ 0.9, which is the value used for the bound plotted in the figure. This not only is consistent with expectations for multivariable models containing CAM noise (32) but also presents a slightly weaker constraint than the single-variable CAM bound plotted in Fig. 2, which has r = 0. The fact that much of the observed data in Fig. 2 satisfy the stronger constraint may therefore be of further significance and deserves to be explored in future work.

Slope of δ¹⁸O versus δ¹³C fluctuations

We estimate the slope of the δ¹⁸O fluctuations versus δ¹³C fluctuations using reduced major axis regression, which is the appropriate choice when both variables contain uncontrolled errors (65). The relationship between δ¹⁸O and δ¹³C changes throughout the Cenozoic has been considered previously (13); however, here, we focus on the sub-million-year fluctuations (i.e., with the long-term trend removed). Figure 4 shows that the model produces δ¹⁸O − δ¹³C slopes, consistent with pre-Pliocene observations. Scatterplots of the observational data, together with the corresponding regression lines, are shown in the Supplementary Materials; it is worth noting that the sign of the slope reverses at the start of the Pliocene, providing further evidence for a switch in the coupling of the climate and the carbon cycle at this time (29).

Adblock test (Why?)

"cycle" - Google News
August 12, 2021 at 01:17AM
https://ift.tt/2VNunAK

Asymmetry of extreme Cenozoic climate–carbon cycle events - Science Advances
"cycle" - Google News
https://ift.tt/32MWqxP
https://ift.tt/3b0YXrX

Bagikan Berita Ini