Stochastic Paleoclimatology: Modeling the EPICA Ice Core Climate Records

We analyze and model the stochastic behavior of paleoclimate time series and assess the implications for the coupling of climate variables during the Pleistocene glacial cycles. We examine 800 kyr of carbon dioxide, methane, nitrous oxide, and temperature proxy data from the EPICA Dome-C ice core, which are characterized by 100~kyr glacial cycles overlain by fluctuations across a wide range of time scales. We quantify this behavior through multifractal time-weighted detrended fluctuation analysis, which distinguishes near red-noise and white-noise behavior below and above the 100~kyr glacial cycle respectively in all records. This allows us to model each time series as a one-dimensional periodic non-autonomous stochastic dynamical system, and assess the stability of physical processes and the fidelity of model-simulated time series. We extend this approach to a four-variable model with linear coupling terms, which we interpret in terms of the interrelationships between the time series. Methane and nitrous oxide are found to have significant destabilizing influences, while carbon dioxide and temperature have smaller stabilizing influences. We draw conclusions about causal relationships in glacial transitions and the climate processes that may have facilitated these couplings, and highlight opportunities to further develop stochastic modeling approaches.


I. INTRODUCTION A. Background & Motivation
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 [1].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 [2].Furthermore, some examples of glacial termination data contradict the hypothesis, on the basis that changes in temperature precede their putative cause of changing insolation [3], and hypothesis testing shows that the 100 ky eccentricity cycle specifically does not significantly influence glacial transitions [4].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][6][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 [8]: One of the most striking features of the 100 ky cycle is its pervasiveness, both geographically and within the various climatic subsystems.It dominates ice mass (and sea level), temperature, carbonate accumulation, upwelling, and carbon dioxide content of the atmosphere.This pervasiveness guarantees that (in the words of Laurent Labeyrie) "everything is correlated with everything", which makes it difficult to deduce mechanisms from proxy records.and by Imbrie et al. [9]: Dozens of explanations have been suggested (section 4).Some models explain the cycle as a free, selfsustaining oscillation with no Milankovitch forcing [e.g., Saltzman and Maasch, 1988].In models of this type, the 100-ky cycle is forced by internal climate system processes so that its phase is arbitrary with respect to eccentricity.Other models explain the cycle as a nonlinear interaction between orbitally forced responses (in the 23-and 41-ky bands) and the internal dynamics of the atmosphere, oceans, ice sheets, and lithosphere [e.g., Maasch and Saltzman, 1990;Gallée et al., 1992].In these, the phase of the 100-ky cycle is orbitally influenced.
For a recent review, the reader is referred to Riechers et al. [10].
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 are 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 [11], 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 [12].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 [e.g., [13][14][15][16].
Paleoclimate analyses have examined causal relationships among paleoclimate data using various approaches, such as comparing prediction quality via convergent cross-mapping [17], quantifying time lag between carbon dioxide and temperature at glacial transitions [18], calculating information flow among variables [19], multivariate autoregressive modeling [20] or using the generalized Fluctuation-Dissipation Relation noted above [11].A multifractal method related to that described here was used by Shao and Ditlevsen [21] 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 [10].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 nonautonomous stochastic dynamical systems theory.In the spirit of other stochastic dynamical systems theory approaches in climate science [22][23][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 [25] using a multifractal analysis method [26] that allows us to identify the color of the noise in the record, and which colors characterize the dynamics over which time scales.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 versus 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 Section 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 Section 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.

Background
We analyze four 800 ky paleoclimate time series extracted from ice cores drilled at Dome C in Antarctica by the EPICA project.We examine records of carbon dioxide (CO 2 ), methane (CH 4 ), nitrous oxide (N 2 O), 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 (δD) 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 [27].
As shown in Figure 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 Figure 2(a).We see the strong 10 −5 yr −1 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 2.5 × 10 −5 and 4.25 × 10 −5 yr −1 , 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 [1].Therefore, we filtered this external forcing by subtracting from each time series a running average (see Section I B 2 for more details).In Figure 2(b) we show the frequency spectrum of each time series after applying the high-pass filter, and we observe that whereas the 10 −5 yr −1 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.

Data preparation
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 [28], particularly in the presence of large gaps in the dataset.Moreover, other studies [29] have confirmed that interpolation generally does not impact 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 interpolating, 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 Section III.

A. Background
We employ multifractal time-weighted detrended fluctuation analysis (MFTWDFA) [26] 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 [30] 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 N  2 for a dataset of length N , as opposed to N 4 in MFDFA.We have used this MFTWDFA in previous work to extract the role of fluctuations in the dynamics of exoplanet detection, sea ice cover and global climate proxy data [31][32][33][34].In order to make this presentation reasonably selfcontained we outline the algorithm presently.

B. Algorithm
To calculate the fluctuation function, we construct a nonstationary profile Y (i) of the original time series X i , as As noted above, in order to work with data evenly spaced in time, we interpolate Y (i) using the modified Akima method.
Next, for each timescale s in the data, the interpolated profile is detrended by removing behavior on timescales longer than that considered.This is done with a point-by-point approximation using weighted linear regression in a window of size s around each point.The weights used in the local linear regression incorporate the intuition that points closer in time are more closely correlated than points farther away in time.Therefore, this continuously weighted fit smoothly captures the local mean and we determine the coefficients for the weighted fit, β, at each point by solving where the elements of the weight matrix W are defined as 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 s.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, 2N s segments are created, where N s = int(N/s) and N is the number of points in the original series X i .
For each timescale s, the variance of the data about the mean is computed up and down the profile, using for ν = [1, N s ], and where ν is the index of the moving time window of size s.Finally, we obtain the fluctuation function, F q (s), as where q denotes the statistical moment.

C. Results: Data Analysis
For q = 2 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, h(2), and the noise type [30].One can relate h(2) to the slope of the power spectrum β as h(2) = (1 + β)/2.For a white noise process β = 0, and hence h(2) = 1/2.For a red-noise process β = 2, and hence h(2) = 3/2.Thus, the varying slopes of the fluctuation function curves demonstrate the different dynamical processes operating on various time scales.
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 a) s < 10 4.6 s > 10  6] 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 (10 3.2 to 10 4.6 years), the fits of the fluctuation functions give h(2) ≈ 3/2, the Hurst exponent of a red-noise process.In the time span between 125 and 400 ky (10 5.1 to 10 5.6 years), we find h(2) ≈ 1/2, 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 the previous section, after which we applied MFTWDFA analysis, which is shown in Figure 3(b).The glacial-scale slopes fall from 1/2 to 0 and the crossover is shifted toward shorter timescales.A Hurst exponent of zero indicates 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 towards 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 towards shorter time scales in the figure.The Hurst exponent for time scales 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 time scales.

A. Background
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 [35].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 time scales 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 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 a(k), b(k), and N (k), 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 (Figure 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 (Figure 3) to find such a peak at a timescale small enough that the time series exhibit red noise dynamics.The largest such frequency is approximately 4.25 × 10 −5 year −1 , or a period of about 23.5 ky, which corresponds to the Milankovitch cycle for the combined effects of axial and apsidal precession [1].

B. One-variable model
We begin with a one-dimensional non-autonomous Ornstein-Uhlenbeck model in order to reproduce the behavior of a single time series, which is based on Moon and Wettlaufer [36], but omits the long-term mean background behavior represented by their f (τ ) term.The non-autonomous model for the time series of the variable η i (t) is The periodic deterministic term a i (t) drives mean-reverting drift, and represents the stability of the process: If a i (t) is negative (positive) the system is stable (unstable) and fluctuations decay (grow) in time.This stability is modulated by the noise amplitude, or noise intensity, given by N i (t), which is also periodic and deterministic.Finally, ξ i (t) is uncorrelated Gaussian white noise, resulting in a red-noise stochastic process, which models the behavior of the paleoclimate data.We solve for a i (t) and N i (t) using a modified version of the procedure described by Moon and Wettlaufer [36] as follows.We consider a time series with M periods and a resolution of T points within each period of length P .Thus, the periodic function a i (t) in Eq. ( 7) is defined as a i (t) = a i ([t/∆t] mod T ), where [.] is the integer part and ∆t = P/T .The same is the case for N i (t).We determine the a i (k) ∀ k ∈ [1, T ] from the analytic solution of the model as where the approximate periodic variance and autocorrelation of the time series of the data X i are defined as respectively.Finally, combining this formulation of a i (k) with the model Langevin equation, the expression for N i (k) is where Step by step details are given in the Supplementary Information of [36].

C. Four-variable model
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 [37] to incorporate CO 2 , CH 4 , N 2 O, 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.
The system of four equations is where η i (t) is the i-th time series, a i (t) is the deterministic stability term, N i (t) is the noise amplitude term, and b ij (k) is the linearized diffusive coupling term representing influence of η j (t) on η i (t), as is common across a wide variety of systems (see Othmer and Scriven [38], Levin [39], Kopell and Howard [40], and Krause et al. [41], for just a few of many examples).
It is significantly simpler to find the a i (k) and b ij (k) in this case, because we can solve four matrix systems for them directly.Each system is constructed by separately multiplying one of the model equations by each of the η i (t), and then taking the ensemble average, resulting in where E xy = ⟨η x (t)η y (t)⟩ and D xy = ⟨ dηx dt η y (t)⟩, which we solve for the coefficients.
Finally, to find the N i (k) we multiply each equation by its corresponding η i (t + ∆t) and take the ensemble average to obtain

D. Results: Model
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.

Deterministic stability coefficients
The deterministic stability in the one-variable model, Eq. ( 7), is controlled by the coefficient a(t), and the deterministic net stability of the four-variable model, Eq. ( 13), is controlled by a i,4D (net) = a i,4D − j b ij .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.

Coupling coefficients and noise amplitude coefficients
The stability and coupling coefficients for the four-variable model are shown in Figures 4(a).The magnitude and sign of the coupling coefficients reflects the interactions between the processes.For simplicity, consider only a two variable system, so that Eq. ( 13) becomes as in Moon and Wettlaufer [37], where the coupling coefficient b 12 (t) (b 21 (t)) represents the influence of variable 2 on variable 1 (variable 1 on variable 2).Therefore, the signs of the coupling coefficients characterize the direction of the influence that a pair of variables have upon each other, and the magnitude characterizes the strength of that interaction.For example, when the coupling coefficient b 12 (t) is positive (negative) then the process represented by variable 2 suppresses (enhances) the growth of variable 1.Thus, in the canonical connotation of stability (instability) viz., the local in time decay (growth) of a variable, a positive (negative) coupling coefficient has a stabilizing (destabilizing) influence on the variables to which it is coupled.Clearly, if b 12 (t) > 0 and b 21 (t) < 0 then variable 2 suppresses the growth of variable 1 and variable 1 enhances the growth of variable 2. [42] Finally, we note that Smale [43] showed that the deterministic form of Eq. ( 16) (i.e., N i (t)ξ i (t) = 0) is a structurally stable global oscillator, that is, apart from a closed set of measure zero, it has a nontrivial periodic attracting solution as t → ∞.The addition of the noise terms N i (t)ξ i (t) 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 Figures 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 exhibit very similar dynamics across all variables.For example, there is one significant peak in the middle of each period, with the exception of the twopeak structure of the one-variable model coefficient for temperature (Figure 4b).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.

FIG. 4: (a)
Coupling coefficients for the four-dimensional model, the one-dimensional model stability and the four-dimensional model net stability, which is defined as a i,4D (net) = a i,4D − j b ij .(b) Noise amplitude coefficients for the one-and four-dimensional models.

Model Interpretation
The principal points at this juncture are as follows.Across all proxies, the time average coupling coefficients for CO 2 and ∆T 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 CH 4 and NO 2 are on average negative, but with a smaller magnitude than those for CO 2 and ∆T.Thus, CH 4 and NO 2 have a weakly destabilizing effect.Finally, the deterministic stability coefficients are all negative.
Clearly, the model captures the canonical strength of the CO 2 and ∆T covariation, and the positive feedback of that covariation on CH 4 and NO 2 .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 [44].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 CH 4 and NO 2 emissions, which increase under a warming climate [45,46], these partially offset the cooling effect of the uptake of carbon by land [47].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 [e.g., 48].We note, however, that this is principally due to their abundances in the atmosphere relative to CO 2 rather than the intrinsic properties of the gases [49].Therefore, the modulation of the CO 2 and ∆T covariation by the warming of the terrestrial biosphere and the associated emission of CH 4 and NO 2 , is consistent with the model presented here.Namely, the relative magnitudes and signs of the coefficients are such that we view the CO 2 and ∆T covariation as the stochastic version of the Smale [43] global oscillator discussed in §III D 2, whose detailed evolution is influenced by the weakly destabilizing dynamics of CH 4 and NO 2 .

Model fidelity
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 Figure 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 Figure 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 Figure 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 towards 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.

Response functions
Knowledge of the model coefficients allows us to construct the linear response matrix function, R(τ ; t), which identifies the causal relations between each time series considered.For the model we study here, R(τ ; t) can be FIG.7: Matrix elements of the time periodic response function obtained from the data plotted as a function of time.
The response function is given by Eq. (17).
written in terms of the time-dependent correlation matrix, C(τ ; t), also called the persistence, according to Baldovin et al. [11] as follows: where the matrix elements of the time-periodic correlation matrix are defined as with M and T defined in Section III B. By expressing the time-dependent correlation function in terms of the model coefficients, we can write the following expression for the time-periodic response function with K ii (t) = a i (t) and K ij (t) = b ij (t).
In Figure 7 we show the temporal behavior of the matrix elements of the response function constructed using the four different data sets, and in Figure 8 we show the matrix elements obtained from the model coefficients.Despite the inherent noise in Figure 7, which results from averaging over only M = 34 points (see Equation 17), we observe qualitative agreement with Figure 8. Specifically, we observe an overall stronger causal relationship from CO 2 to ∆T than from ∆T to CO 2 , consistent with the findings of Baldovin et al. [11].However, the strength of these causal links varies throughout the period, with certain time windows displaying a stronger causal relationship from CO 2 to ∆T, while others exhibit the reverse relationship.Furthermore, we observe that CH 4 and NO 2 have a negligible influence on ∆T and CO 2 , while ∆T and CO 2 have a strong causal link to, and hence strong influence on, CH 4 and NO 2 .Importantly, this analysis of the response functions is consistent with the model interpretation discussed in §III D 3.

IV. CONCLUSION
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 intervariable 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 time scale 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 [50] 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 §III D 3. Stocker et al. [47] note that the contemporary terrestrial biosphere mitigates anthropogenic climate change by acting as a carbon sink, which compensates approximately 30% 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 versus 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 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 [51] and Persson [52] used multivariate Granger causality tests in their analyses of the EPICA ice core data to show that CO 2 , ∆T and CH 4 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 §III D and mentioned throughout is our treatment of only four variables, which may themselves be coupled to others.For example, the analysis of CH 4 in the EPICA ice core data by Loulergue et al. [53] indicates that the connection between ice-sheet volume and Antarctic temperature and CH 4 millennial variability, and the relationship proposed in the literature between them, fails to capture millennial CH 4 events in the early glacial phases.Thus, the coupling between CH 4 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 CO 2 and ∆T, 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][46][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.

FIG. 1 :
FIG. 1: Carbon dioxide, methane, temperature, and nitrous oxide time series from the EPICA ice core record.(a) Original time series, (b) normalized time series of the fluctuations relative to the slowly-varying mean.

FIG. 2 :
FIG. 2: Power spectra of EPICA time series for (a) original data and (b) data fluctuations relative to the slowly-varying mean behavior.
4.5 for the shorter-timescale regime and s > 5.1 for the longer-timescale regime.(a) Original EPICA time series, (b) forcing-converted time series, based on linear regression fits to the fluctuation function below and above the 100 ky glacial cycle crossover.timescales over which we can examine the noise behavior.By looking at Figure 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, s = [3.2,4.6] for the original dataset and s = [3.2,3.6] for the fluctuations, and s = [5.1, 5.

FIG. 3 :
FIG. 3: Logarithmic plots of the fluctuation functions for (a) original EPICA time series and (b) time series of fluctuations, with the slowly-varying behavior removed.Coarse wide gray lines show regressions fitted to approximately straight segments of the fluctuation function that correspond to distinct regimes of colored noise behavior, and vertical dotted grey lines show the 23 ky periodicity used in the later modeling section.

FIG. 5 :
FIG. 5: Comparison (from -600 ky to -400 ky) of (a) forcing data time series with simulated time series generated from (b) one-variable and (c) four-variable models.

FIG. 6 :
FIG. 6: Comparison of (a) periodic standard deviations, (b) probability density functions, and (c) autocorrelation functions between the data and one-variable (top rows) and four-variable (bottom rows) models.From left to right the columns are CO 2 , CH 4 , ∆T and NO 2 .

FIG. 8 :
FIG. 8: Matrix elements of the time dependent response function obtained from the model coefficients plotted as a function of time.The response function is given by Eq. (19).

TABLE I :
Scaling exponent estimates from MFTWDFA fluctuation function slopes for the two colored-noise regimes, using s < 10