Stochastic paleoclimatology: Modeling the EPICA ice core climate records

The Earth’s Quaternary glacial cycles are characterized by noisy processes with differing dynamics across timescales, but similar large-scale periodic behavior among different climate variables that correspond to the glacial cycles of the Pleistocene. An area of particular interest in paleoclimate dynamics is the origin of the 100 ky cycle. The canonical explanation for glacial cycle pacing is the Milankovitch hypothesis, which attributes it to periodic changes in Earth’s orbital parameters.¹ Namely, because variations in Earth’s eccentricity, obliquity, and precession change the distance and angle of incident insolation to the planet’s surface over time, the resulting temperature changes are thought to drive the variations in greenhouse gas concentrations that are seen during glacial cycles. The Milankovitch cycle for eccentricity has an approximately 100 ky period, matching the glacial cycle periodicity.

However, this hypothesis is the subject of great scrutiny, as evidence for it is typically based on pattern matching between the insolation and paleoclimate datasets. It is unclear why glacial cycles would be paced by eccentricity because it is the weakest of the Milankovitch cycles, as the hypothesis itself does not explain what kinds of mechanisms could amplify this small signal into one that dominates glacial pacing.² Furthermore, some examples of glacial termination data contradict the hypothesis, on the basis that changes in temperature precede their putative cause of changing insolation,³ and hypothesis testing shows that the 100 ky eccentricity cycle specifically does not significantly influence glacial transitions.⁴ Indeed, a substantial challenge involves clearly identifying the internal climate mechanisms and feedbacks governing glacial cycles and in particular, the interactions between paleoclimate variables. The community understands that many of the physical and chemical mechanisms that can facilitate these interactions, including the greenhouse effect, ocean carbon uptake, carbon rock weathering, soil nitrogen release, and permafrost melt,^5–7 and seeks understanding of which processes may have dominated paleoclimate dynamics and hence may underlie the pace of glaciations. The issues were succinctly summarized by Berger and Wefer,⁸ 

For a recent review, the reader is referred to Riechers et al.¹⁰ 

Although our goal here is not to put forth a new theory for glacial pacing, we are interested in understanding the stochastic dynamics, noise characteristics, and causal relationships among several key paleoclimate proxies that accompany glaciations. To that end, the development of models that reproduce multiscale stochastic dynamics and elucidate causal interactions among climate processes is our focus. Many common statistical methods, for example, computing the covariance, can tell us the strength of the relationship between two variables, but cannot reveal the direction of cause and effect within that relationship, nor whether one process stabilizes or destabilizes another. This problem can be addressed using a generalized Fluctuation–Dissipation Relation,¹¹ which is able to identify causal links between the processes, but cannot reveal stabilizing and destabilizing relationships between them. Global climate models simulate interactions in the climate system by numerically integrating conservation laws throughout the atmosphere and ocean and incorporating the influence of forcings and parameterization of relatively small-scale processes.¹² However, they often cannot reproduce the variability, small-scale structure, and long duration typical of climate time series due to the limited treatment of, and intermodel differences between, subgrid-scale processes as well as the processing power needed to run such models over long time periods.^13–16 

Paleoclimate analyses have examined causal relationships among paleoclimate data using various approaches, such as comparing prediction quality via convergent cross-mapping,¹⁷ quantifying time lag between carbon dioxide and temperature at glacial transitions,¹⁸ calculating information flow among variables,¹⁹ multivariate autoregressive modeling,²⁰ or using the generalized Fluctuation–Dissipation Relation noted above.¹¹ A multifractal method related to that described here was used by Shao and Ditlevsen²¹ to study the different scaling properties of interglacial and glacial climates using a wide range of data, including Antarctic and Greenland ice cores. They found the Holocene record to be monofractal, and the glacial record to be multifractal, and concluded that the glacial climate has a longer persistence time and stronger nonlinearities.

These approaches reach a variety of conclusions about the dominant causal direction among temperature and greenhouse gases and about the validity of the Milankovitch hypothesis.¹⁰ This motivates new approaches of examining causal paleoclimate relationships.

Our approach here is to use a stochastic data analysis and modeling method involving colored noise and non-autonomous stochastic dynamical systems theory. In the spirit of other stochastic dynamical systems theory approaches in climate science,^22–24 we can characterize the random variability of climate processes that are not captured by deterministic models.

We first quantify the types of noise present in the time series²⁵ using a multifractal analysis method²⁶ that allows us to identify the color of the noise in the record, and which colors characterize the dynamics over which timescales. This is essential for climate time series that exhibit both significant noise behavior and timescale separation so that we can assess how the dynamics differ on shorter vs longer timescales.

We then model paleoclimate time series as Ornstein–Uhlenbeck processes, consisting of periodic, nonautonomous Langevin equations that treat both the deterministic behavior and stochastic variability of the record. We apply these models to carbon dioxide, methane, nitrous oxide, and temperature proxy time series and assess their performance by computing the response function for each pair of variables.

Finally, we interpret the physical significance of the resulting model coefficients and the response functions and pursue their implications for the interactions among these paleoclimate variables. With this in hand, we can simulate the time series and assess their fidelity relative to the original data through a variety of statistical metrics.

The structure of this paper is as follows. In Sec. II, we apply multifractal time-weighted detrended fluctuation analysis (MFTWDFA) to paleoclimate ice core records and quantify the nature of the fluctuations found therein. We introduce and apply our Ornstein–Uhlenbeck models in Sec. III and examine their properties and fidelity through statistical comparisons and computation of the response functions. Having extracted the noise types present in these paleoclimate records, we reproduce the behavior and examine the causal relationships between proxies using simple stochastic models. We conclude with a discussion of the implications of these analyses for the last 800 000 years of Earth’s climate history.

We analyze four 800 ky paleoclimate time series extracted from ice cores drilled at Dome C in Antarctica by the European Project for Ice Coring in Antarctica (EPICA) project. We examine records of carbon dioxide (⁠⁠), methane (⁠⁠), nitrous oxide (⁠⁠), and a proxy for temperature. Greenhouse gas concentrations are estimated by direct measurement of air bubbles trapped in the ice, while temperature is reconstructed from the deuterium proxy (⁠⁠) and presented as the change in temperature compared to the 1950 average global temperature. Direct measurements of nitrous oxide are supplemented by measurements of nitrous oxide artifacts where direct measurements were not possible. The four datasets use the EDC3 chronology, based on snow accumulation, flow modeling, and independent age markers, to estimate the correspondence between core depth and age.²⁷ 

As shown in Fig. 1, the 100 ky periodic glacial cycles are clearly observed in the time series. Additionally, however, the data also exhibit a complex, noisy structure across timescales. We first examine this structure by computing the frequency spectrum of all four time series, as shown in Fig. 2(a). We see the strong peak corresponding to the period of 100 ky, the same period as the eccentricity cycle (lower-frequency peaks are not reliable because of the lack of data.), and, as noted above, is a feature that is the focus of a great deal of debate and research.

Furthermore, we observe other relevant peaks near and ⁠, corresponding approximately to the 41 ky obliquity cycle and 23.5 ky net precession cycle resulting from the combination of axial and apsidal precession. These peaks are typically attributed to the presence of external astronomical forcing in all the time series, which makes them highly correlated and, consequently, makes their causality relationships extremely difficult to unravel.¹ Therefore, we filtered this external forcing by subtracting from each time series a running average (see Sec. I B 2 for more details). In Fig. 2(b), we show the frequency spectrum of each time series after applying the high-pass filter, and we observe that whereas the peak is significantly reduced, the other two are increased. Clearly, we have filtered (not expunged) the time varying external forcing and hence there remains an associated footprint in the time series, which must be taken into account. This motivates building a non-autonomous stochastic system to model the filtered data and to examine causal relationships.

In order to use the approach described above, we interpolate the time series to an evenly spaced temporal resolution. We interpolate to match as closely as possible the lowest-resolution dataset—nitrous oxide, with 912 points—while also splitting the time domain into the 34 equal periods, resulting in 25 points per period and a spacing of approximately 929 years between points. This constant spacing does mean that in the original time series, multiple points may be interpolated into some time gaps, but we confirmed that this is relatively uncommon and most time gaps in the original data are on the scale of this interpolation gap—the main exception being a large time gap in the nitrous oxide time series, which is an unavoidable limitation of the EPICA dataset. For interpolation, we utilize the Akima method, which eliminates unrealistic overshoots introduced by other interpolation methods, such as the cubic spline,²⁸ particularly in the presence of large gaps in the dataset. Moreover, other studies²⁹ have confirmed that interpolation generally does not impact the results of statistical analysis. Due to the 20 ky gap in the nitrous oxide record from 260 to 240 ky, an artificial data point was added to the nitrous oxide time series at 250 ky using linear interpolation in that domain to better constrain Akima interpolation for our analysis.

After interpolation, we removed the slow-varying mean behavior, as we focus on modeling the smaller scale fluctuations. We applied a Gaussian filter with a smoothing filter using a characteristic time-window three times the time increment used for interpolation. This approach resulted in an optimal filtering of slow fluctuations, while maintaining the fast fluctuations. We subtracted this mean behavior from the interpolated time series to obtain fluctuations around the mean. Subsequently, we normalized the fluctuation time series so that each has a standard deviation of unity and, thus, can be modeled comparably. The distributions of the resulting fluctuation time series are nearly Gaussian, which supports our modeling approach described in Sec. III.

We employ multifractal time-weighted detrended fluctuation analysis (MFTWDFA)²⁶ to extract the scaling dynamics and fluctuation structure in the EPICA paleoclimate time series. This method quantifies the fluctuations in the time series around the mean behavior across timescales present in the data through the fluctuation function defined below. The approach enables us to draw conclusions about the dominant statistical fluctuations as a function of timescale. If the fluctuations in a time series are colored noise, the fluctuation function will scale exponentially over increasingly large timescales, and the particular value of the scaling exponent, referred to as the Hurst exponent, corresponds to the color. Therefore, a log–log plot of the fluctuation function is a straight line over the range of time in which the data exhibit a particular colored fluctuation behavior, and the slope of this line will be the corresponding Hurst exponent. A power spectrum analysis could accomplish the same goal as MFTWDFA of quantifying colored noise, but the multifractal approach provides a clearer and more accurate description of the complex multiscale nature of the paleoclimate data and, in particular, the crossover times between regimes of noise behavior.

MFTWDFA builds on other detrended fluctuation analysis methods such as MFDFA³⁰ by introducing a smoother computation of the mean behavior of the data on each timescale. In MFDFA, a piecewise polynomial fit to the profile of the data is used. In MFTWDFA, a time-weighted linear regression in a moving window provides a continuous estimate of the mean behavior at each timescale, leading to a fluctuation function that shows crossover times between noise regimes more clearly. Furthermore, MFTWDFA allows us to extract information about the nature of fluctuations at timescales up to for a dataset of length ⁠, as opposed to in MFDFA.

We have used this MFTWDFA in the previous work to extract the role of fluctuations in the dynamics of exoplanet detection, sea ice cover and global climate proxy data.^31–34 In order to make this presentation reasonably self-contained, we outline the algorithm presently.

We then start at the beginning of the profile and split the data into intervals with an equal number of points, whose total time range corresponds to the timescale ⁠. Accounting for the possibility that a portion of the profile remains, the same operation is reversed beginning at the end of the profile. In this manner, segments are created, where and is the number of points in the original series ⁠.

For ⁠, we fit straight line segments to the logarithmic plots of the fluctuation functions, which allows us to determine the Hurst exponents of the time series at different timescales. We use the second moment due to the simplicity of the correspondence between the Hurst exponent, and the noise type.³⁰ One can relate to the slope of the power spectrum as ⁠. For a white noise process, ⁠, and hence ⁠. For a red-noise process, ⁠, and hence ⁠. Thus, the varying slopes of the fluctuation function curves demonstrate the different dynamical processes operating on various timescales.

Ideally, to demonstrate robust scaling behavior, the straight-line slope segments should span as many orders of magnitude as possible. However, the length and resolution of the time series set upper and lower bounds for the timescales over which we can examine the noise behavior.

By looking at Fig. 3(a), we can see that all four of these climate variables are governed by similar stochastic dynamics below and above the glacial cycle timescale of 100 ky. Furthermore, the fluctuation functions for the four original EPICA time series clearly show two distinct regimes of colored noise behavior. We fit straight lines to these two regions, for the original dataset and for the fluctuations, and for the long-timescale side, and find their slopes in order to quantify their noise types. In the time span between 1.5 and 40 ky (⁠ to years), the fits of the fluctuation functions give ⁠, the Hurst exponent of a red-noise process. In the time span between 125 and 400 ky (⁠ to years), we find ⁠, exhibiting white-noise behavior. The fluctuation function structure for nitrous oxide differs slightly from the others at smaller timescales, with some subtler crossovers rather than the single slope seen in the other datasets. However, for the purpose of our modeling, we are mainly concerned with the fluctuation function slopes for the data after the high-pass filter has removed the slowly varying behavior.

We then normalized the data and applied the high-pass filter as described in Sec. I B 2, after which we applied MFTWDFA, which is shown in Fig. 3(b). The glacial-scale slopes fall from to 0 and the crossover is shifted toward shorter timescales. See Table I for a summary. A Hurst exponent of zero indicates the lack of scaling behavior on longer scales, and thus mean-reversion behavior of the time series. This can be attributed to the fact that the filtered time series describes a stochastic process that decays toward a constant position rather than a slowly time-varying signal. As the length of the smoothing windows is reduced, higher frequencies are removed from the data, leading to the observed shift in the transition point toward shorter timescales in the figure. The Hurst exponent for timescales below the length of the applied smoothing window of approximately 3 ky remained unaltered since the high-pass filter does not remove frequencies below the smoothing window. As a result, we can confidently use a multidimensional non-autonomous Ornstein–Uhlenbeck process to model the data on these timescales.

Climate time series can be modeled via simple stochastic processes if there is a clear separation between short and long timescales with distinct dynamics, and if the short-term processes can be modeled as random walks.³⁵ Such a modeling approach is of interest because it incorporates the small-scale random fluctuations typical of climate processes into a modeling framework that can be run over much longer timescales than can be achieved by current global climate models. We are also interested in the simplicity of such models, which allows us to determine the parameters of interest, such as stability and noise amplitude, analytically, and we can easily introduce coupling functions that can illuminate the nature of the interactions in the climate system.

Our results from MFTWDFA justify the use of such a framework to model and analyze the EPICA paleoclimate time series. The difference in short-term and long-term dynamics for all four variables shows a clear separation of timescales between the sub-glacial and super-glacial periods. Moreover, even after filtering the long-term glacial cycle behavior, the short-term sub-glacial behavior is a nonstationary, approximately red noise, time series.

The stochastic model we employ extends the Ornstein–Uhlenbeck process to a non-autonomous periodic system with a separation of timescales. An Ornstein–Uhlenbeck process is the overdamped limit of the Langevin equation describing the Brownian motion, when the particle experiences the restoring influence of a local quadratic potential. Thus, the potential causes the dynamics to be mean-reverting. We can add to this canonical stationary process a longer-timescale forcing term to represent the slowly varying mean behavior of the glacial cycles. Therefore, such a model can appropriately represent the way our paleoclimate time series fluctuates around this slowly varying mean behavior.

Our model coefficients—the drift and the noise amplitude terms—are time-dependent and periodic. We seek to model the strong Milankovitch frequencies present in the paleoclimate time series, and this periodicity allows us to derive the model coefficients from the periodic statistics of the data. We choose the period of the model coefficients ⁠, ⁠, and ⁠, described below, based on the power spectra of our climate variables, Milankovitch cycle periodicities, and the timescale-separated noise structure revealed in MFTWDFA. We compare the power spectra (Fig. 2) of our four time series to identify a frequency peak, corresponding to a Milankovitch period, that is common across all four datasets. We also examine the MFTWDFA fluctuation functions (Fig. 3) to find such a peak at a timescale small enough that the time series exhibits red noise dynamics. The largest such frequency is approximately  year⁠, or a period of about 23.5 ky, which corresponds to the Milankovitch cycle for the combined effects of axial and apsidal precession.¹ 

We extend this one-dimensional model in order to treat multiple time series together, introducing coupling terms that represent the influence that each time series has on the others. These coupling terms allow us to make first-order estimates of the primary direction of influence between time series variables, information beyond what measures like the covariance can provide. We develop a four-variable model based on Moon and Wettlaufer³⁷ to incorporate ⁠, ⁠, ⁠, and temperature time series and their couplings to each other. We note that this type of system can be extended to an arbitrary number of variables.

We apply this modeling approach to the EPICA paleoclimate time series to derive and interpret the stability, coupling, and noise coefficients for each of the four variables in the coupled system. We quantify the model fidelity by using these coefficients to simulate artificial time series and compare them with the original time-series data.

The deterministic stability in the one-variable model, Eq. (7), is controlled by the coefficient ⁠, and the deterministic net stability of the four-variable model, Eq. (13), is controlled by (net) ⁠. We note that these coefficients are comparable across models in each variable, as we would expect since they are treating the same process. However, the four-variable model is slightly more negative across all four processes, showing that the couplings between the processes enhance the overall stability. In both models, methane and nitrous oxide are more stable than are carbon dioxide and temperature, and hence their deterministic drift drives them more strongly toward the long-term mean behavior.

Finally, we note that Smale⁴³ showed that the deterministic form of Eq. (17) [i.e., ] is a structurally stable global oscillator, that is, apart from a closed set of measure zero, it has a nontrivial periodic attracting solution as ⁠. The addition of the noise terms in our model simply “smears out” the attracting solution to a degree that depends on the noise amplitude. We return to this below.

We see in Fig. 4(a) that, in the main, the coupling coefficients connecting carbon dioxide and temperature to each other or to methane or nitrous oxide have positive signs, and hence act to suppress the growth of these variables. On the other hand, the coupling coefficients of methane or nitrous oxide are negative throughout, indicating that they act to enhance the growth of other variables. The generally larger magnitudes of the coupling coefficients for carbon dioxide and temperature indicate their dominant control. Of course, the overall dynamics depends on all of the terms in Eq. (13).

The periodic behavior of the noise terms of the one- and four-variable models exhibits very similar dynamics across all variables. For example, there is one significant peak in the middle of each period, with the exception of the two-peak structure of the one-variable model coefficient for temperature [Fig. 4(b)]. However, across all variables, the one-variable noise amplitude is consistently larger than the four-variable value. This is simply because the coupling terms in the latter provide additional sources of fluctuations, and hence each variable’s individual noise amplitude compensates by contributing a smaller amount of noise, thereby maintaining the same overall noise level between models.

The principal points at this juncture are as follows. Across all proxies, the time average coupling coefficients for and are positive, and hence their mutual interactions are stabilizing. However, as seen in Fig. 4(a), depending on time, one can be larger than the other in an approximately periodic manner so that the mutual stabilization is time dependent. In contrast, the coupling coefficients for and are on average negative, but with a smaller magnitude than those for and ⁠. Thus, and have a weakly destabilizing effect. Finally, the deterministic stability coefficients are all negative.

Clearly, the model captures the canonical strength of the and covariation, and the positive feedback of that covariation on and ⁠. Contemporary studies show that, in general, warming-induced methane-climate feedbacks are positive, with the principal contributors being atmospheric methane lifetime and biogenic emissions from wetlands and permafrost.⁴⁴ Such feedbacks are complicated by the fact that the terrestrial biosphere presently acts as a partial compensatory carbon sink of global emissions. Indeed, because the terrestrial biosphere is responsible for substantial fractions of and emissions, which increase under a warming climate,^45,46 these partially offset the cooling effect of the uptake of carbon by land.⁴⁷ Moreover, because adding nitrous oxide (methane) is about 200 (20) times more effective at increasing global temperatures as adding equal amounts of carbon dioxide, small fluctuations in the emissions of nitrous oxide and methane could be amplified into large effects on climate.⁴⁸ We note, however, that this is principally due to their abundances in the atmosphere relative to rather than the intrinsic properties of the gases.⁴⁹ Therefore, the modulation of the and covariation by the warming of the terrestrial biosphere and the associated emission of and ⁠, is consistent with the model presented here. Namely, the relative magnitudes and signs of the coefficients are such that we view the and covariation as the stochastic version of the Smale⁴³ global oscillator discussed in Sec. III D 2, whose detailed evolution is influenced by the weakly destabilizing dynamics of and ⁠.

We use the model coefficients computed from the EPICA time series and run our one- and four-variable models forward in time using a standard Euler method. This generates artificial time series with statistics and noise behavior that should match those of the original detrended time series. Figure 5 shows that both models reproduce the general appearance of the four time series. Next, we compare key statistical metrics to quantify how well our simulations reproduce different aspects of the proxy data.

In Fig. 6(a), we compare the periodic standard deviations of the data and models for each variable and find that the four-variable model is superior to the one-variable model in that it reproduces the overall magnitude and periodic shape of the standard deviation quite well for all of the time series. The probability distribution functions also compare favorably, although we observe that the model fits a Gaussian distribution to the slightly non-Gaussian observations, as shown in Fig. 6(b). Thus, while both models reproduce the observational mean and standard deviation (within 3% of the observational statistics in all cases), they do not reproduce the skewness and the kurtosis.

In Fig. 6(c), we compare the autocorrelation functions, which are less well reproduced than the other statistics. Whereas the one- and four-variable models both capture some of the oscillations in the autocorrelation function of the data, neither model reproduces the magnitude of the negative minimum of the data, nor the decay rate toward that minimum value. Here again, apart from some model approximations, which may not capture the full complexity of nonlinear processes in these paleoclimate time series, there may be many additional variables in the observed system that couple to those four in the observed record, but cannot be reflected in the four we treat in the model. However, we note that the four-variable model reproduces the rate of decay in autocorrelation and some of the negative values better than does the one-variable model, indicating an important role of the coupling coefficients. Nonetheless, we view this behavior of the autocorrelation as a weakness in the predictive power of our approach.

In Fig. 7, we show the temporal behavior of the matrix elements of the response function constructed using the four different data sets, and in Fig. 8 we show the matrix elements obtained from the model coefficients. Despite the inherent noise in Fig. 7, which results from averaging over only points [see Eq. (18)], we observe qualitative agreement with Fig. 8. Specifically, we observe an overall stronger causal relationship from to than from to ⁠, consistent with the findings of Baldovin et al.¹¹ However, the strength of these causal links varies throughout the period, with certain time windows displaying a stronger causal relationship from to ⁠, while others exhibit the reverse relationship. Furthermore, we observe that and have a negligible influence on and ⁠, while and have a strong causal link to, and hence a strong influence on, and ⁠. Importantly, this analysis of the response functions is consistent with the model interpretation discussed in Sec. III D 3.

The Earth’s paleoclimate underwent periodic but noisy 100 ky cycles of glaciation and deglaciation over the last 800 ky, which are clearly visible in time series data for carbon dioxide, methane, nitrous oxide, and temperature obtained from the EPICA ice core. We used a multifractal method to study these time series and extract the types of colored noise that characterize them across scales, as well as the times at which there is a crossover between behaviors, in a more precise way than the usual spectral slope analysis allows. This allowed us to adopt and extend previous non-autonomous stochastic models to represent each paleoclimate time series individually, and then as a coupled system, taking into account the time-dependent structure of their deterministic and stochastic dynamics.

Our combined approach produces observationally consistent simple stochastic dynamical models. We extracted the timescale-separated colored noise regimes in the data, and computed and interpreted the stability, noise, and inter-variable couplings through non-autonomous Ornstein–Uhlenbeck models. These coupling coefficients demonstrate the directionality and magnitude of the stabilizing effects of interactions between these climate variables, providing insight into the multiple timescale dynamics of the climate.

A central finding of our stochastic treatment is that carbon dioxide and temperature have stabilizing influences on each other and on methane and nitrous oxide, but the latter two have a weakly destabilizing influence on each other and on carbon dioxide and temperature. The strong co-variation between carbon dioxide and temperature has long been the signature of glacial cycles, but with the perennial question regarding which variable drives the other (see, e.g., Cuffey and Vimeux⁵⁰ and references therein). Both the stochastic model coefficients and the response functions show this carbon dioxide and temperature “pulse” of the climate system, but with a time-dependence of which one has a controlling influence. The weakly destabilizing influence of methane and nitrous oxide is due to the positive feedback—enhanced emissions—of the terrestrial biosphere to warming as discussed in Sec. III D 3. Stocker et al.⁴⁷ note that the contemporary terrestrial biosphere mitigates anthropogenic climate change by acting as a carbon sink, which compensates approximately of global carbon dioxide emissions. Moreover, given the efficacy of methane and nitrous oxide as greenhouse gases, and the destabilizing influence we have identified our approach, it is clear that the carbon dioxide and temperature pulsing of glaciations is modulated by the terrestrial biosphere. Keeping in mind that we have only modeled four proxies, we note that the asymmetry between stadials and interstadials, with the abrupt warming vs more gradual cooling, is consistent with our analysis. The high latitude terrestrial biosphere is snow and ice covered during a stadial and the ice-albedo feedback exhibits hysteresis. Thus, abrupt ice loss is accompanied by the abrupt release of methane and nitrous oxide and thereby facilitates rapid warming. During the warm interstadial, slow terrestrial carbon uptake facilitates cooling until sufficient snow and ice cover suppresses terrestrial emissions driving the climate into a stadial.

On the one hand, Kang and Larsson⁵¹ and Persson⁵² used multivariate Granger causality tests in their analyses of the EPICA ice core data to show that ⁠, ⁠, and all “Granger cause” each other. Namely, their analysis strongly rejects the null hypothesis that any of these three variables does not cause the other. On the other hand, one of the important caveats discussed in Sec. III D and mentioned throughout is our treatment of only four variables, which may themselves be coupled to others. For example, the analysis of in the EPICA ice core data by Loulergue et al.⁵³ indicates that the connection between ice-sheet volume and Antarctic temperature and millennial variability, and the relationship proposed in the literature between them, fails to capture millennial events in the early glacial phases. Thus, the coupling between and other climate variables not treated here is non-trivial and not simply reflected in the coupling coefficients. Therefore, although the signature of glacial cycles is generally principally associated with the covariation of and ⁠, our results suggest that the interactive-coupled role of other greenhouse gases is important in the timing of these cycles. This realization is of course not new,^45–47 but the point here is that it is cast in a framework that is much simpler to use than a comprehensive climate model.

The approach described here constitutes a modest step in quantifying paleoclimate glacial dynamics using simple stochastic modeling techniques. Natural advances in the modeling framework would include, among others, nonlinear coupling between variables, nonlinear multiplicative, and/or correlated noise. However, having examined only four paleoclimate observables, which may be coupled to many other variables, a clear next step is to introduce additional variables into our modeling framework. Clearly, this requires incorporation of more proxy variables, thereby increasing the complexity of the coupled model, but such an approach is nonetheless vastly simpler than using comprehensive global climate models. Finally, because our approach reproduces key statistical dynamical quantities, it can in principle act as a constraint for comprehensive global climate models across a range of observationally accessible epochs.

N.D.B.K. and J.S.W. gratefully acknowledge support from the Yale University. L.T.G. and J.S.W. gratefully acknowledge support from the Swedish Research Council (Vetenskapsrådet) under Grant No. 638-2013-9243. Nordita is partially supported by Nordforsk.

The authors have no conflicts to disclose.

N. D. B. Keyes: Investigation (equal). L. T. Giorgini: Investigation (equal). J. S. Wettlaufer: Conceptualization (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.