Presence or absence of stabilizing Earth system feedbacks on different time scales

Stabilizing feedbacks, in principle, should affect how the typical amplitude of fluctuations within a system changes with time scale (21). To show how this would work, we take a purposely simplified perspective of the Earth system in which the only variable of interest is globally averaged surface temperature, T. This simplification is appropriate for a first attempt at extracting information about long-term Earth system feedbacks directly from data of past fluctuations; furthermore, as we will show, it is already sufficient for obtaining useful insight.

Two simple “end-member” scenarios for this simplified view are displayed in Fig. 1. Scenario A is the classic established model of climate variability in the absence of stabilizing feedbacks: a random walk (22–24). This assumes that slowly evolving components of the Earth system retain an aggregate “memory” of the fast-evolving components that accumulates approximately randomly (22). In that case, temperature evolution would be described by the following stochastic differential equation where η(t) is a Gaussian white noise forcing and a is a constant. In this model, the root mean square temperature fluctuation ΔT_rms occurring on a time scale Δt is proportional to Δt^1/2 (equivalent to red noise; see Materials and Methods). Many climate time series exhibit this scaling behavior (22–26), and the ability to reproduce it is part of the model’s appeal. Throughout this paper, we will often refer to the scaling exponent (1/2 in this case) as H.

Scenario B not only is the same as scenario A but also includes a stabilizing feedback with characteristic (i.e., e-folding) time scale τ [also known as an Ornstein-Uhlenbeck process (21)]

On time scales Δt ≪ τ, the feedback term is negligible and the root mean square fluctuation still scales as Δt^1/2. However, the feedback damps correlations for time scales Δt ≫ τ, and the root mean square fluctuation then scales as Δt^−1/2 (Materials and Methods). Further, aggregating multiple stabilizing feedback processes on different time scales can yield apparent power laws ΔT_rms ∝ Δt^H for any −1/2 < H < 1/2 (27–29) (see also Materials and Methods).

The real Earth system is of course much more complicated than this. There are a vast range of processes on a vast range of time scales that are not explicitly accounted for. Nevertheless, as the pioneering work by Hasselmann (22) showed, in complex systems such as Earth’s climate, the combined effects of many deterministic processes can be aggregated by the slower components of the system to yield statistics essentially like a random walk (scenario A above). Thus, the η(t) in Eqs. 1 and 2 can be considered to already account for many of these processes; the explicit feedback term in Eq. 2 just means that there is a dominant stabilizing feedback on a time scale τ.

Long-term feedbacks in the real Earth system do not necessarily act directly on temperature. For example, of the two mentioned in the Introduction, the silicate weathering feedback responds directly to temperature and the carbonate compensation feedback does not. Nevertheless, if long-term temperature variability is driven at least in part by variability in atmospheric CO₂, any feedback that helps stabilize CO₂ is indirectly helping to stabilize temperature.

A final point needs to be made regarding the possibility of periodic forcings and resonances. On geologic time scales, climate is forced by periodic oscillations in Earth’s orbital parameters (30, 31); these forcings, if powerful enough, could be expected to create a peak in fluctuation amplitudes similar to that in scenario B (Fig. 1). The same would be true if the Earth system had an intrinsic tendency to oscillate at a certain time scale. A case study for both would be Plio-Pleistocene glacial variability, and this will be worth addressing once we take a look at the data.

We calculate the root mean square temperature fluctuation ΔT_rms as a function of time scale Δt for five different paleotemperature time series (Materials and Methods). We consider four benthic foraminiferal δ¹⁸O records (32–36) and one compilation of isotopic temperatures from Antarctic ice cores (37): Between them, they resolve fluctuations on time scales spanning more than five orders of magnitude. Specifically, “fluctuations” are defined using Haar wavelets (20, 38). Considering a time series of temperature, T(t), the Haar fluctuation ΔT over a time interval Δt is defined as the difference between the average values of the time series over the first and second halves of the interval; this is described schematically in Fig. 2 and discussed further in Materials and Methods. We use it because it is simple, accurately measures scaling behavior (38), and is straightforwardly applied to unevenly sampled paleoclimate time series (20). It also highlights the physically important difference between fluctuations growing with scale (H > 0) or shrinking with scale (H < 0).

The results of our analysis are shown in Fig. 3; some power-law scalings (with fixed exponents H) are added as guides for interpretation. A previous analysis by Lovejoy (20) suggested the existence of three regimes that are relevant here: a “climate” regime on time scales below about 80 ka in which fluctuations increase with time scale, a “macroclimate” regime in which fluctuations decrease with time scale, and a “megaclimate” regime above about 500 ka in which fluctuations increase with time scale again. Our analysis paints a similar picture but with some key differences.

On time scales shorter than about 4 ka and longer than about 400 ka, fluctuations increase with time scale: H ≃ 0.5, similar to a random walk and consistent with scenario A. Between 4 and 400 ka, the behavior depends on what interval the data cover. Datasets that contain exclusively Plio-Pleistocene variability (i.e. the last) show a clear peak at a few tens of thousands of years and a strongly decreasing regime beyond this; this forms the basis of the regime classification by Lovejoy (20) noted above. However, our analysis reveals that throughout the rest of the Cenozoic, these fluctuations consistently obeyed H ≃ 0—that is, their amplitude is essentially time scale independent. The anomalous Plio-Pleistocene peak and the regime with rapidly decreasing fluctuation amplitudes beyond it likely record the rapid periodic transitions between glacial and interglacial states, rather than evidence regarding stabilizing feedbacks (see Materials and Methods for a further discussion).

Following the previous section and Fig. 1, the fact that H is much less than 0.5 in this intermediate regime strongly suggests that stabilizing feedbacks have exerted dominant control over Earth’s surface temperature on time scales between 4 and 400 ka. We emphasize how remarkable it is that the amplitude of the typical root mean square fluctuation in global temperature is essentially constant across two orders of magnitude in time scale. While our analysis cannot conclusively show which feedbacks were responsible, we can make inferences by comparing the time scales to those of various known or hypothesized feedbacks: This is what we will do in the Discussion. To aid this, Fig. 3 also shows the approximate time scales of important Earth system feedbacks in this regime, as well as their likely signs (see Materials and Methods for details).

To make clear how multiple feedbacks in a complex system can create a regime with time scale–independent ΔT_rms as in Fig. 3, and to help develop a more specific interpretation of the three regimes shown in the data, we expand on the stochastic models discussed earlier. Specifically, we consider Earth’s surface temperature T to be the sum of multiple stochastic processes, some with stabilizing feedbacks (e.g., scenario B) and some without (scenario A). Mathematically, we let where $\stackrel{\dot{}}{f_i}=-f_i/{\mathrm{\tau }}_i+a_i{\mathrm{\eta }}_i(t)$, $\stackrel{\dot{}}{r}(t)=a_n{\mathrm{\eta }}_n$, and η_i are independent Gaussian white noise forcings (discussed further in Materials and Methods). Last, a_n < a_i for all i < n, meaning that variability due to the random walk r(t) grows more slowly than that of the other processes. A key property of this model is that the stabilizing feedbacks have only partial control—in other words, they only stabilize part of the system, and there can still be undamped variability at other scales. The real Earth system shares this property: if it did not, paleoclimate records would exhibit no variability at all on long time scales.

As an example, we choose partial stabilizing feedbacks on time scales of 1, 10, and 100 ka (τ₁, τ₂, and τ₃, respectively), numerically simulate Eq. 3 for 200 Ma, and analyze fluctuations using the same algorithm that we applied to the real data. Results are shown in Fig. 4; the general behavior of the observations is well reproduced. On short time scales (<τ₁), fluctuations grow like a random walk with H ≃ 0.5 and then have essentially time scale–independent amplitudes in the regime in which the feedbacks are active. On long time scales (>τ₃), the undamped stochastic variability (reflecting the partial nature of the feedbacks) takes over, and fluctuations again grow like a random walk. Theory predicts that this kind of behavior occurs for a wide range of possible models and parameter values (Materials and Methods): In all cases, the position of the intermediate regime is determined by the range of time scales of stabilizing feedbacks.

The identification of the anomalous 4- to 400-ka regime is a novel confirmation that stabilizing feedbacks with characteristic time scales in this regime have been a dominant control on Earth’s surface temperature. To understand which mechanisms were likely responsible, we can compare this time scale range to the previously proposed time scales for different stabilizing feedbacks.

Of immediate interest is the consistency of this regime with the ∼100-ka time scale proposed for the silicate weathering feedback (7, 8). We suggest that this is strong observational evidence for the importance of silicate weathering as a climate stabilizer. Through this, it further supports the widely used steady-state assumption (4), existing models of the long-term effects of anthropogenic CO₂ emissions (6, 8), and the idea that the weathering feedback should play a key role in planetary habitability (5).

The fact that the nongrowing regime seems to start at time scales as small as 4 ka suggests that other shorter–time scale stabilizing feedbacks were also important. One obvious candidate is ocean mixing: The ocean can help damp temperature fluctuations due to its large thermal inertia, and full equilibration is achieved on a time scale of a few thousand years (39). Another possibility on a ∼10-ka time scale is CaCO₃ equilibration (6, 8, 9), which could indirectly stabilize temperature through its effect on atmospheric CO₂. Other feedbacks potentially active at this time scale include vegetation and land ice (see Fig. 3); however, these are likely both destabilizing (mathematically positive) feedbacks (40) and hence would not have been responsible for stabilization.

What is the origin of the increasing fluctuation amplitude beyond 400 ka? Following the theory explained above, the random walk–like growth (H ≃ 0.5) should mean that there are no dominant stabilizing feedbacks in the system on these time scales. Yet, if current thinking is at all accurate, the silicate weathering feedback should still be active on these time scales: It is not inherently time scale limited. What then is going on?

One possible resolution is the following. Consider Earth’s “weathering curve” (41), interpreted here as the dependence of the silicate weathering flux, F_si, on Earth’s surface temperature. Neglecting changes in organic carbon oxidation or burial, a steady state is established when F_si is equal to the volcanic flux F_volc of carbon into the surface environment. Because the weathering curve has a positive slope (weathering increases with temperature), we obtain the familiar stabilizing feedback that tends to drive the system toward a steady state.

Nevertheless, the weathering curve itself may change over time (41) due to changes either in the surface carbon cycle’s physical attributes [such as the amount and properties of exposed weatherable rock (13, 42–44)] or in the mechanisms constituting the feedback itself [e.g., due to biological evolution in land plants (1)] (see also Fig. 3B). This will lead to slow “quasistatic” changes in the surface temperature (45), even while the carbon cycle remains in steady state with respect to input and output fluxes. We suggest that it is precisely this class of changes that lead to fluctuations increasing again at the longest time scales.

Figure 5 summarizes this schematically. Imagine that the weathering curve moves upward, for example, due to an increase in weatherability; then, the new steady state will move to a lower surface temperature. On time scales of hundreds of thousands of years, the weathering feedback will damp fluctuations toward the steady state. Yet, on longer time scales, weathering will act as a forcing, and the steady state itself will move. The H ≃ 0.5 scaling beyond 400 ka then suggests that the steady state moves in an undamped way. In other words, while silicate weathering is a stabilizing feedback bringing the system to a steady state, there are no stabilizing feedbacks on the million-year–time scale motion of that steady state itself.

This is of course a highly simplified picture of weathering. We have ignored the effects of changes in organic carbon oxidation and burial, and are considering factors such as CO₂, topography, vegetation types, precipitation, and rock types only implicitly (by arguing that they change the weathering curve). Nevertheless, we suggest that the basic reasoning regarding a weathering-established steady state that moves in an undamped way is likely independent of these details. This could and should be tested using a more detailed carbon cycle model.

Last, there is one other possibility that deserves mention: that the increasing fluctuation amplitudes at the longest time scales are due to other destabilizing feedbacks. While there are no obvious candidate mechanisms for these feedbacks on multi-million year time scales, the data at present cannot rule this out. This would also be very interesting to pursue further.

The fact that global temperature fluctuations continue to grow like a random walk at the longest time scales has major implications for understanding the long-term habitability of Earth and other Earth-like planets. There is a long-standing debate (14–16) over the extent to which Earth’s observed billion-year habitability is a product of stabilization (for example, due to the weathering feedback) or a product of chance. The predominant view has been that the weathering feedback is responsible for this long-term habitability, and this stabilization is a key part of the habitable zone concept used to search for life on other planets (5, 46).

We have shown that the observations are inconsistent with a dominant stabilizing feedback on the longest time scales and suggested that those fluctuations arise because of weathering acting as a climate forcing (for example, when tectonic processes change the availability of weatherable rocks). Another option is that fluctuations grow on long time scales because of unknown destabilizing feedbacks. In either case, the key question is: Are there any mechanisms in the Earth system that prevent these kinds of fluctuations from eventually driving surface temperature into an uninhabitable regime? If there are none, it would follow that chance may have played a nonnegligible role in Earth’s continued habitability, and that other Earth-like planets with an active carbonate-silicate cycle and in the conventional habitable zone may not necessarily be as accommodating to life over long periods of time as has previously been expected. Obtaining and analyzing well-calibrated, higher-resolution paleotemperature records spanning longer stretches of geologic time, as well as improving our understanding of tectonic evolution and its climatic consequences on time scales of many millions of years (47), should provide further insights.

It has long been recognized that climate time series on various time scales exhibit self-similar “scaling” behavior (20, 22, 24, 25, 48, 49). A process x(t) is considered to exhibit self-similarity if where H is the self-similarity exponent and the equality is in terms of probability distributions. For these processes, the power spectrum S(ω) ∝ ω^−β, where β ≃ 2H + 1 (50). When β ≃ 1 (i.e., H ≃ 0), this is the widely studied “1/f noise” (51). Observed climate time series often exhibit well-defined time scale regimes in which β and H take on different values (20, 49).

To begin to understand the physical origin of this scaling, a simple null model without feedbacks considers climate fluctuations as a random walk (22). In the limit of infinitesimal step sizes, this is the Wiener process (21), which has probability distribution

Relating this to Eq. 4, one can show that this process is self-similar with H = 1/2 (i.e., β = 2). In contrast, white noise, which is the long-time limit of Eq. 2, has β = 0 (21) (i.e., H = − 1/2).

Superposing multiple processes with stabilizing feedbacks at different time scales can create apparent scaling exponents in the range −1/2 < H < 1/2 (27). For the benefit of the reader, the Supplementary Text demonstrates this explicitly. In the literature on 1/f noise, H = 0 has also been associated with a continuous log-uniform distribution of feedback time scales (28, 29).

Self-similarity in real data can be measured either in real space or in frequency space. In real space, this can be done through a fluctuation function Δx(Δt) [“structure function” in the field of turbulence (52)]. Ideally, this function would obey

One possible choice of fluctuation function is the simple difference but this only accurately reflects scaling behavior (i.e., behaves according to Eq. 7) in the regime 0 < H < 1 (38). In this work, we define Δx(Δt) as the Haar fluctuation: the difference between the time series averaged over the first and second halves of the interval Δt (38, 53). This accurately reflects scaling behavior in the range −1 < H < 1 (38) and is additionally desirable because of its conceptual and computational simplicity. In particular, it is straightforward to measure the scaling behavior for unevenly sampled time series without interpolation. The details of the algorithm we use are described in the “Algorithm for calculating fluctuation time scale dependence” section.

We analyze four benthic δ¹⁸O datasets: the Cenozoic compilations first introduced by Zachos et al. (32) and updated by Zachos et al. (33), the Cenozoic composite record of Westerhold et al. (36), the orbitally tuned Plio-Pleistocene compilation of Lisiecki and Raymo (34), and the non-orbitally tuned Pleistocene compilation of Huybers (35). Temperature values are inferred from δ¹⁸O following the calibrations proposed by Hansen et al. (54) and slightly redefined by Westerhold et al. (36), which take into account changes in the ice volume contribution to δ¹⁸O and the relationship between deep-ocean and surface temperatures throughout the Cenozoic. This conversion is based on absolute δ¹⁸O data, and thus the Huybers (35) time series, which has long-term averages removed, has to have a constant offset added to match the others (although the specific value does not affect our results). Deep-ocean temperature T_do, following the calibrations discussed above, is given by

Surface temperature T is then calculated from the deep-ocean temperature according to

While the above conversion is used for the sake of accuracy, it is worth noting that the qualitative results remain the same even if a single linear conversion constant relating δ¹⁸O and T is used for all datasets.

Last, we also include in our analysis the ice core temperature compilation of Parrenin et al. (37). We divide the fluctuations by a factor of 2 to approximately account for high-latitude amplification (20), although again it is important to emphasize that the details of this correction do not affect our conclusions.

We calculate the root mean square fluctuation for each paleotemperature time series, ΔT_rms, using an interpolation-free algorithm based on that proposed by Lovejoy (20): All of our codes are made freely available at (55) and https://github.com/arnscheidt/stabilizing-earth-system-feedbacks. For an unevenly sampled time series described by two vectors t = [t₀, t₁, t₂…t_n] and T = [T₀, T₁, T₂…T_n], we define the Haar fluctuation at position j of size k: i.e., the difference between the average value of the second k/2 and the first k/2 data points. This can be implemented efficiently using cumulative sums.

We consider all even k ranging from 0 to n. Although we could calculate all possible ΔT(j, k) for the data, ΔT(j₁, k) is negligibly different from ΔT(j₂, k) as long as k ≫ ∣ j₁ − j₂∣. Therefore, for each k, we calculate ΔT(j, k) at intervals of ak, where we choose a = 0.5: This gives the algorithm n log n time complexity instead of n². Because the interval k is divided in two in terms of indices, but not in terms of time elapsed, we discard any Δx(j, k) with $\mathrm{ϵ}<\frac{t_{j+k/2}-t_j}{t_{j+k}-t_j}<1-\mathrm{ϵ}$ for some ϵ: The choice of ϵ defines a balance between robustly allowing for unevenly spaced data but ensuring that the extracted ΔT(Δt) remain meaningful. Following Lovejoy (20), we use ϵ = 0.25.

Last, we calculate Δt = t_{j + k} − t_k for each ΔT(j, k) and average the ΔT(Δt)² over evenly spaced bins in log space (four bins per order of magnitude). The data points in Fig. 3 are the bin centers, and the root mean square fluctuation ΔT_rms is given by taking the square root of the averaged ΔT(Δt)². At the extremes (very small or large Δt), there begin to be much fewer data points; we therefore truncate the data where the number of data points per bin are a factor b smaller than the maximum (we use b = 5).

For the purposes of Fig. 3, we have additionally truncated some of the shortest–time scale fluctuations (Δt < 4 ka) from the data of Zachos et al. (33) and Westerhold et al. (36). They are anomalously large compared to those in the more recent higher-resolution datasets, as well as the results when the 20-Ma segments of the Westerhold et al. (36) data are studied individually. The latter do not consider the time interval from 5 Ma to present, suggesting that this effect arises only due to the data in that interval. Meanwhile, the higher-resolution datasets from this same period [Lisiecki and Raymo (34), Huybers (35), and Parrenin et al. (37)] represent averages over multiple records, while the Zachos et al. (33) and Westerhold et al. (36) data do not. Therefore, we suggest that the anomaly probably reflects contributions to the δ¹⁸O signal from sources other than global temperature (regional-scale variability, diagenesis, etc.), which are not of interest here. The 20-Ma segments shown in Fig. 3 have not been additionally truncated in this manner, so none of this affects our conclusions.

Figure 3 shows that the datasets spanning only the Plio-Pleistocene (∼5 Ma to present) (34, 35, 37) behave somewhat differently in the intermediate regime than the datasets spanning the entire Cenozoic (32, 36): There is a more pronounced peak at time scales of a few tens of thousands of years, and a regime in which fluctuations decrease very rapidly. What is the origin of this difference? A likely solution is that it simply reflects the onset of the Plio-Pleistocene glacial cycles, which feature marked transitions between different climate states on time scales of tens of thousands of years. This variability would produce an anomalous peak in the averaged fluctuation amplitudes near this time scale and thus correspondingly steeper slopes on either sides of the peak.

Repeated periodic transitions such as the glacial cycles are a confounding effect when attempting to observe signatures of stabilizing feedbacks in the data. Nevertheless, we know that this kind of behavior is limited to the Plio-Pleistocene: There should be no such confounding effect in other time periods. As shown in Fig. 3, the data from the rest of the Cenozoic show the same consistent pattern of fluctuations not growing between time scales of about 4 to 400 ka, providing strong evidence of dominant control by stabilizing feedbacks.

The long-term Earth system feedbacks and their time scales shown in Fig. 3 are loosely taken after Rohling et al. (56), who made distinctions between annual, decadal, century, millennial, multimillennial, and million year time scales. We have also included the signs of the feedbacks (i.e., positive/destabilizing or negative/stabilizing) where this is clear [following, e.g., (4, 6, 7, 40, 56)]. The time scales are still only intended as approximate; here, we offer some additional justification for the more specific values shown. We emphasize that the upper time scale limit for a given feedback does not mean that the underlying process is not operating on longer time scales: It simply means that the feedback has reached steady state.

The upper limit of the land ice feedback can be taken approximately as the time scale on which Plio-Pleistocene deglacations occur (10⁴ years). More specifically, we can take it as the time scale at which we see the peak in the Plio-Pleistocene temperature fluctuations (i.e., 3 ×10⁴ years). The upper limit of the vegetation feedback is taken approximately as 10³ years, based on the observation of strong climate-vegetation correlations on millennial time scales (57). The upper limit of the stabilizing feedback due to ocean mixing is taken as a few thousand years, based on simulations showing that substantial tracer disequilibrium can persist for at least 2000 years (39). Last, we split carbon cycle feedbacks into three categories: those on time scales shorter than 5 ka (which may be stabilizing or destabilizing), CaCO₃ equilibration from 5 to 50 ka (6, 8), and silicate weathering operating between 50 ka and a few million years (4, 7, 8).

To make concrete the relationship between stabilizing feedbacks and the time scale dependence of fluctuations as shown in Fig. 3, we consider stochastic models of long-term temperature variability. Stochastic models have long been used in the study of glacial cycles (24) but have only recently begun to be applied to deep time climate problems (58–60).

Considering Eq. 3, for n = 2 (i.e., one stabilizing feedback), one can show that (Supplementary Text)

This already begins to qualitatively reproduce the three regimes seen in the data. For n > 2 (i.e., with more stabilizing feedback processes), things become more complicated; yet, it is not hard to obtain behavior like

The first crossover time scale τ₁ is that of the fastest feedback, and the slow crossover time scale τ_s is typically determined by the feedback with the slowest time scale together with the amplitude of the slow random walk r(t). Depending on the choice of parameters, it is now possible to reproduce any behavior of the kind seen in Fig. 3.

For the example in Fig. 4, n = 4, τ₁ = 1 ka, τ₂ = 10 ka, τ₃ = 100 ka, a₁ = 0.03 K year^−1/2, a₂ = 0.0085 K year^−1/2, a₃ = 0.0027 K year^−1/2, a₄ = 0.0015 K year^−1/2. As can be seen in the figure, the intermediate regime has H ≃ 0. The 200-Ma numerical simulation was carried out using an Euler-Maruyama algorithm and the Julia package DifferentialEquations.jl (61).

We emphasize that while we have chosen the noise sources η_i to be uncorrelated for simplicity, our results should, in principle, be independent of this. For example, the result that multiple stochastic feedback processes can be superimposed to create scaling regimes with −1/2 < H < 1/2 holds both for uncorrelated noise sources (Supplementary Text) and for correlated ones (27).