# Stochastic Paleoclimatology: Modeling the EPICA Ice Core Climate Records

N. D. B. Keyes<EMAIL_ADDRESS>Program in Applied Mathematics, Yale
University, New Haven, Connecticut 06520, USA  Department of Earth & Planetary
Sciences, Yale University, New Haven, Connecticut 06520, USA  L. T. Giorgini
<EMAIL_ADDRESS>Nordic Institute for Theoretical Physics, Royal
Institute of Technology
and Stockholm University, Stockholm 10691, Sweden  J.S. Wettlaufer
<EMAIL_ADDRESS>Program in Applied Mathematics, Yale University, New
Haven, Connecticut 06520, USA  Department of Earth & Planetary Sciences, Yale
University, New Haven, Connecticut 06520, USA  Nordic Institute for
Theoretical Physics, Royal Institute of Technology
and Stockholm University, Stockholm 10691, Sweden  Department of Physics, Yale
University, New Haven, Connecticut 06520, USA

###### Abstract

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 kiloyears of carbon dioxide,
methane, nitrous oxide and temperature proxy data from the EPICA Dome-C ice
core, which are characterized by 100 ky 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
ky glacial cycle respectively in all records. This allows us to model each
time series as a one-dimensional periodic nonautonomous 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 intervariable coupling terms, which we interpret in terms of possible
interrelationships among the four time series. Within the framework of our
coupling coefficients, we find that carbon dioxide and temperature act to
stabilize each other and methane and nitrous oxide, whereas the latter two
destabilize each other and carbon dioxide and temperature. We also compute the
response function for each pair of variables to assess the model performance
by comparison to the data and confirm the model predictions regarding
stability amongst variables. Taken together, our results are consistent with
glacial pacing dominated by carbon dioxide and temperature that is modulated
by terrestrial biosphere feedbacks associated with methane and nitrous oxide
emissions.

> Paleoclimate time series of greenhouse gases and temperature show Earth’s
> periodic but noisy glacial transitions over the last 800,000 years, which
> are widely attributed to periodic changes in orbital forcing, but are still
> not well understood. Here, we apply a multifractal analysis method to four
> time series from the EPICA ice core to understand its colored-noise
> structure across timescales, showing that they are characterized by near-red
> noise on subglacial timescales and near-white noise on glacial timescales.
> Informed by this result, we model the time series as stochastic processes,
> first individually and then as a linearly coupled system, to extract
> stability and noise coefficients and assess interactions among the
> variables. We examine the model coupling coefficients and compute the
> response function for each pair of variables, to reveal stabilizing and
> destabilizing relationships of varying strengths between them, suggesting
> potential causal relationships in climate transitions.

##  I Introduction

###  I.1 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, self-sustaining 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 non-autonomous 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.

Figure 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. Figure 2:
Power spectra of EPICA time series for (a) original data and (b) data
fluctuations relative to the slowly-varying mean behavior.

###  I.2 The EPICA dataset

####  I.2.1 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 ($\text{CO}_{2}$), methane ($\text{CH}_{4}$), nitrous oxide
($\text{N}_{2}\text{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 ($\delta
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}$ $\text{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\times 10^{-5}$ and
$4.25\times 10^{-5}$ $\text{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.2.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}$ $\text{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.

####  I.2.2 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.

##  II Multifractal Time-Weighted Detrended Fluctuation Analysis

###  II.1 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 $\frac{N}{2}$ for a dataset of
length $N$, as opposed to $\frac{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 self-
contained we outline the algorithm presently.

###  II.2 Algorithm

To calculate the fluctuation function, we construct a nonstationary profile
$Y(i)$ of the original time series $X_{i}$, as

$Y(i)\equiv\sum_{k=1}^{i}\Big{(}X_{k}-\bar{X}\Big{)},\hskip
14.22636pti=1,...,N.$ (1)

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, $\hat{\beta}$, at each point by solving

$(X^{T}WX)\hat{\beta}=X^{T}Wy,$ (2)

where the elements of the weight matrix $W$ are defined as

$w_{ij}=\begin{cases}\Big{(}1-(\frac{i-j}{s})^{2}\Big{)}^{2},&|i-j|\leq s\\\
0,&\text{otherwise.}\end{cases}$ (3)

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}=\text{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

$V(\nu,s)=\frac{1}{s}\sum_{i=1}^{s}\Big{[}Y([\nu-1]s+i)-\hat{y}([\nu-1]s+i)\Big{]}^{2},$
(4)

for $\nu=[1,N_{s}]$, and

$V(\nu,s)=\frac{1}{s}\sum_{i=1}^{s}\Big{[}Y(N-[\nu-N_{s}]s+i)-\hat{y}(N-[\nu-
N_{s}]s+i)\Big{]}^{2},$ (5)

for $\nu=[N_{s}+1,2N_{s}]$, where $\nu$ is the index of the moving time window
of size $s$.

Finally, we obtain the fluctuation function, $F_{q}(s)$, as

$F_{q}(s)=\Big{[}\frac{1}{2N_{s}}\sum_{\nu=1}^{2N_{s}}\\{V(\nu,s)\\}^{\frac{q}{2}}\Big{]}^{\frac{1}{q}},$
(6)

where $q$ denotes the statistical moment.

###  II.3 Results: Data Analysis

Figure 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.

a)  $s<10^{4.6}$ $s>10^{5.1}$ CO2 1.50 0.46 CH4 1.30 0.19 $\Delta\mathrm{T}$
1.48 0.18 N2O 1.19 0.49 b)  $s<10^{3.6}$ $s>10^{5.1}$ CO2 1.47 0.01 CH4 1.26
0.01 $\Delta\mathrm{T}$ 1.30 0.01 N2O 1.34 0.00

Table 1: Scaling exponent estimates from MFTWDFA fluctuation function slopes
for the two colored-noise regimes, using $s<10^{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.

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
$\beta$ as $h(2)=(1+\beta)/2$. For a white noise process $\beta=0$, and hence
$h(2)=1/2$. For a red-noise process $\beta=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
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.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)\approx 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)\approx 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.

##  III Stochastic models

###  III.1 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\times 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].

###  III.2 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(\tau)$ term. The non-autonomous model for the time
series of the variable $\eta_{i}(t)$ is

$\frac{d\eta_{i}}{dt}=a_{i}(t)\eta_{i}(t)+N_{i}(t)\xi_{i}(t).$ (7)

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, $\xi_{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/\Delta t]\,\textrm{mod}\,T)$, where $[.]$ is the integer
part and $\Delta t=P/T$. The same is the case for $N_{i}(t)$. We determine the
$a_{i}(k)\leavevmode\nobreak\ \forall\leavevmode\nobreak\ k\in[1,T]$ from the
analytic solution of the model as

$a_{i}(k)\approx-\frac{1}{\Delta t}\frac{S_{i}(k)-A_{i}(k)}{S_{i}(k)},$ (8)

where the approximate periodic variance and autocorrelation of the time series
of the data $X_{i}$ are defined as

$\displaystyle S_{i}(k)$
$\displaystyle\equiv\frac{1}{M-1}\sum_{j=1}^{M}X^{jT+k}_{i}X^{jT+k}_{i}\approx\langle(\eta_{i}(k))^{2}\rangle,\qquad\textrm{and}$
(9) $\displaystyle A_{i}(k)$
$\displaystyle\equiv\frac{1}{M-1}\sum_{j=1}^{M}X^{jT+k}_{i}X^{jT+k+1}_{i}\approx\langle\eta_{i}(k)\eta_{i}(k+1)\rangle,$
(10)

respectively.

Finally, combining this formulation of $a_{i}(k)$ with the model Langevin
equation, the expression for $N_{i}(k)$ is

${N_{i}(k)}=\sqrt{\frac{\langle y_{i}^{2}(k)\rangle}{\Delta t}},$ (11)

where

$y_{i}(k)\equiv\eta_{i}(k+1)-\eta_{i}(k)-a_{i}(k)\eta_{i}(k)\Delta t.$ (12)

Step by step details are given in the Supplementary Information of [36].

###  III.3 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 $\text{CO}_{2}$, $\text{CH}_{4}$, $\text{N}_{2}\text{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

$\begin{split}\frac{\mathrm{d}\eta_{i}(t)}{\mathrm{d}t}&=a_{i}(t)\eta_{i}(t)+N_{i}(t)\xi_{i}(t)+\sum_{j\neq
i}b_{ij}(t)\Big{[}\eta_{j}(t)-\eta_{i}(t)\Big{]}\end{split}$ (13)

where $\eta_{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 $\eta_{j}(t)$ on
$\eta_{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 $\eta_{i}(t)$, and then taking the ensemble average, resulting in

[Eii, Eij\- Eii, Eik\- Eii, Eil\- EiiEij, Ejj\- Eij, Ejk\- Eij, Ejl\- EijEik,
Ejk\- Eik, Ekk\- Eik, Ekl\- EikEil, Ejl\- Eil, Ekl\- Eil, Ell\- Eil] ×[ai(t)
bij(t) bik(t) bil(t) ] = [DiiDijDikDil] ,

(14)

where $E_{xy}=\langle\eta_{x}(t)\eta_{y}(t)\rangle$ and
$D_{xy}=\langle\frac{d\eta_{x}}{dt}\eta_{y}(t)\rangle$, which we solve for the
coefficients.

Finally, to find the $N_{i}(k)$ we multiply each equation by its corresponding
$\eta_{i}(t+\Delta t)$ and take the ensemble average to obtain

$\begin{split}N_{i}^{2}(t)&=\Big{\langle}\eta_{i}(t+\Delta
t)\frac{d\eta_{i}}{dt}\Big{\rangle}-a_{i}(t)\Big{\langle}\eta_{i}(t)\eta_{i}(t+\Delta
t)\Big{\rangle}\\\ &-b_{ij}(t)\Big{[}\langle\eta_{i}(t+\Delta
t)\eta_{j}(t)\rangle-\langle\eta_{i}(t+\Delta t)\eta_{i}(t)\rangle\Big{]}\\\
&-b_{ik}(t)\Big{[}\langle\eta_{i}(t+\Delta
t)\eta_{k}(t)\rangle-\langle\eta_{i}(t+\Delta
t)\eta_{i}(t)\rangle\Big{]}.\end{split}$ (15)

###  III.4 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.

####  III.4.1 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}-\sum_{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.

####  III.4.2 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

$\displaystyle\frac{\mathrm{d}\eta_{1}(t)}{\mathrm{\leavevmode\nobreak\ d}t}$
$\displaystyle=a_{1}(t)\eta_{1}(t)+N_{1}(t)\xi_{1}(t)+b_{12}(t)\left[\eta_{2}(t)-\eta_{1}(t)\right]\qquad\text{and}$
$\displaystyle\frac{\mathrm{d}\eta_{2}(t)}{\mathrm{\leavevmode\nobreak\ d}t}$
$\displaystyle=a_{2}(t)\eta_{2}(t)+N_{2}(t)\xi_{2}(t)+b_{21}(t)\left[\eta_{1}(t)-\eta_{2}(t)\right],$
(16)

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.111We note, however, that a different connotation of stability can
be used [37]. Namely, in the context of the longevity of a climate variable,
when $b_{12}(t)$ is negative (positive), and thereby provides a weak positive
(negative) forcing to variable 1, we can say that variable 2 stabilizes
(destabilizes) the presence in the climate system of variable 1.

Finally, we note that Smale [43] showed that the deterministic form of Eq.
(16) (i.e., $N_{i}(t)\xi_{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\rightarrow\infty$. The addition
of the noise terms $N_{i}(t)\xi_{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 two-
peak 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.

Figure 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}-\sum_{j}b_{ij}$. (b) Noise
amplitude coefficients for the one- and four-dimensional models. Figure 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.

Figure 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 $\mathrm{CO}_{2}$, $\mathrm{CH}_{4}$, $\Delta\mathrm{T}$ and
$\mathrm{NO}_{2}$.

####  III.4.3 Model Interpretation

The principal points at this juncture are as follows. Across all proxies, the
time average coupling coefficients for $\mathrm{CO}_{2}$ and
$\Delta\mathrm{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
$\mathrm{CH}_{4}$ and $\mathrm{NO}_{2}$ are on average negative, but with a
smaller magnitude than those for $\mathrm{CO}_{2}$ and $\Delta\mathrm{T}$.
Thus, $\mathrm{CH}_{4}$ and $\mathrm{NO}_{2}$ have a weakly destabilizing
effect. Finally, the deterministic stability coefficients are all negative.

Clearly, the model captures the canonical strength of the $\mathrm{CO}_{2}$
and $\Delta\mathrm{T}$ covariation, and the positive feedback of that
covariation on $\mathrm{CH}_{4}$ and $\mathrm{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
$\mathrm{CH}_{4}$ and $\mathrm{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 $\mathrm{CO}_{2}$ rather
than the intrinsic properties of the gases [49]. Therefore, the modulation of
the $\mathrm{CO}_{2}$ and $\Delta\mathrm{T}$ covariation by the warming of the
terrestrial biosphere and the associated emission of $\mathrm{CH}_{4}$ and
$\mathrm{NO}_{2}$, is consistent with the model presented here. Namely, the
relative magnitudes and signs of the coefficients are such that we view the
$\mathrm{CO}_{2}$ and $\Delta\mathrm{T}$ covariation as the stochastic version
of the Smale [43] global oscillator discussed in §III.4.2, whose detailed
evolution is influenced by the weakly destabilizing dynamics of
$\mathrm{CH}_{4}$ and $\mathrm{NO}_{2}$.

####  III.4.4 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.

####  III.4.5 Response functions

Knowledge of the model coefficients allows us to construct the linear response
matrix function, $R(\tau;t)$, which identifies the causal relations between
each time series considered. For the model we study here, $R(\tau;t)$ can be
written in terms of the time-dependent correlation matrix, $C(\tau;t)$, also
called the persistence, according to Baldovin  _et al._ [11] as follows:

$R(\tau;t)=C(\tau;t)C^{-1}(\tau;0),$ (17)

where the matrix elements of the time-periodic correlation matrix are defined
as

$C_{ij}(\tau;t)=\frac{1}{M}\sum_{n=1}^{M}\eta_{i}(\tau+t+nT\Delta
t)\eta_{j}(\tau+nT\Delta t),$ (18)

with $M$ and $T$ defined in Section III.2. 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

$R(\tau;t)=\textrm{exp}\left[\int_{\tau}^{\tau+t}{\bf
K}(t^{\prime})\mathrm{d}t^{\prime}\right],$ (19)

with $K_{ii}(t)=a_{i}(t)$ and $K_{ij}(t)=b_{ij}(t)$.

Figure 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). Figure 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).

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
$\mathrm{CO}_{2}$ to $\Delta\mathrm{T}$ than from $\Delta\mathrm{T}$ to
$\mathrm{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
$\mathrm{CO}_{2}$ to $\Delta\mathrm{T}$, while others exhibit the reverse
relationship. Furthermore, we observe that $\mathrm{CH}_{4}$ and
$\mathrm{NO}_{2}$ have a negligible influence on $\Delta\mathrm{T}$ and
$\mathrm{CO}_{2}$, while $\Delta\mathrm{T}$ and $\mathrm{CO}_{2}$ have a
strong causal link to, and hence strong influence on, $\mathrm{CH}_{4}$ and
$\mathrm{NO}_{2}$. Importantly, this analysis of the response functions is
consistent with the model interpretation discussed in §III.4.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 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 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.4.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 $\mathrm{CO}_{2}$, $\Delta\mathrm{T}$ and $\mathrm{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.4 and mentioned throughout is
our treatment of only four variables, which may themselves be coupled to
others. For example, the analysis of $\mathrm{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 $\mathrm{CH}_{4}$ millennial
variability, and the relationship proposed in the literature between them,
fails to capture millennial $\mathrm{CH}_{4}$ events in the early glacial
phases. Thus, the coupling between $\mathrm{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 $\mathrm{CO}_{2}$ and
$\Delta\mathrm{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.

###### Acknowledgements.

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

## Author Declarations

### Conflict of Interest

The authors have no conflicts to disclose.

## References

  * Hays, Imbrie, and Shackleton [1976] J. D. Hays, J. Imbrie, and N. J. Shackleton, “Variations in the Earth’s Orbit: Pacemaker of the Ice Ages,” Science 194, 1121–1132 (1976).
  * Wunsch [2004] C. Wunsch, “Quantitative estimate of the Milankovitch-forced contribution to observed Quaternary climate change,” Quat. Sci. Rev. 23, 1001–1012 (2004).
  * Winograd  _et al._ [1992] I. J. Winograd, T. B. Coplen, J. M. Landwehr, A. C. Riggs, K. R. Ludwig, B. J. Szabo, P. T. Kolesar, and K. M. Revesz, “Continuous 500,000-Year Climate Record from Vein Calcite in Devils Hole, Nevada,” Science 258, 255–260 (1992).
  * Huybers and Wunsch [2005] P. Huybers and C. Wunsch, “Obliquity pacing of the late Pleistocene glacial terminations,” Nature 434, 491–494 (2005).
  * Dean  _et al._ [2018] J. F. Dean, J. J. Middelburg, T. Röckmann, R. Aerts, L. G. Blauw, M. Egger, M. S. M. Jetten, A. E. E. de Jong, O. H. Meisel, O. Rasigraf, C. P. Slomp, M. H. in’t Zandt, and A. J. Dolman, “Methane Feedbacks to the Global Climate System in a Warmer World,” Reviews of Geophysics 56, 207–250 (2018).
  * Williams, Katavouta, and Goodwin [2019] R. G. Williams, A. Katavouta, and P. Goodwin, “Carbon-Cycle Feedbacks Operating in the Climate System,” Current Climate Change Reports 5, 282–295 (2019).
  * Xu-Ri  _et al._ [2012] Xu-Ri, I. C. Prentice, R. Spahni, and H. S. Niu, “Modelling terrestrial nitrous oxide emissions and implications for climate feedback,” New Phytologist 196, 472–488 (2012).
  * Berger and Wefer [2003] W. H. Berger and G. Wefer, “On the Dynamics of the Ice Ages: Stage-11 Paradox, Mid-Brunhes Climate Shift, and 100-ky Cycle,” in  _Earth’s Climate and Orbital Eccentricity: the Marine Isotope Stage 11 Question_ , Geophys. Monogr 137, edited by A. Droxler, R. Z. Poore, and L. H. Burckle (AGU, 2003) pp. 41–59.
  * Imbrie  _et al._ [1993] J. Imbrie, A. Berger, E. A. Boyle, S. C. Clemens, A. Duffy, W. R. Howard, G. Kukla, J. Kutzbach, D. G. Martinson, A. McIntyre, A. C. Mix, B. Molfino, J. J. Morley, L. C. Peterson, N. G. Pisias, W. L. Prell, M. E. Raymo, N. Y. Shackleton, and J. R. Toggweiler, “On the structure and origin of major glaciation cycles. 2. the 100,000-year cycle,” Paleoceanography 8, 699–735 (1993).
  * Riechers  _et al._ [2022] K. Riechers, T. Mitsui, N. Boers, and M. Ghil, “Orbital insolation variations, intrinsic climate variability, and quaternary glaciations,” Clim. Past 18, 863–893 (2022).
  * Baldovin  _et al._ [2022] M. Baldovin, F. Cecconi, A. Provenzale, and A. Vulpiani, “Extracting causation from millennial-scale climate fluctuations in the last 800 kyr,” Scientific Reports 12, 15320 (2022), number: 1 Publisher: Nature Publishing Group.
  * Hansen  _et al._ [1983] J. Hansen, G. Russell, D. Rind, P. Stone, A. Lacis, S. Lebedeff, R. Ruedy, and L. Travis, “Efficient Three-Dimensional Global Models for Climate Studies: Models I and II,” Monthly Weather Review 111, 609–662 (1983).
  * Stone and Risbey [1990] P. H. Stone and J. S. Risbey, “On the limitations of general circulation climate models,” Geophysical Research Letters 17, 2173–2176 (1990).
  * Lopez-Gomez  _et al._ [2020] I. Lopez-Gomez, Y. Cohen, J. He, A. Jaruga, and T. Schneider, “A generalized mixing length closure for eddy-diffusivity mass-flux schemes of turbulence and convection,” J. Adv. Model. Earth Syst. 12 (2020), 10.1029/2020MS002161.
  * Ma  _et al._ [2022] P.-L. Ma, B. E. Harrop, V. E. Larson, R. B. Neale, A. Gettelman, H. Morrison, H. Wang, K. Zhang, S. A. Klein, M. D. Zelinka, Y. Zhang, Y. Qian, J.-H. Yoon, C. R. Jones, M. Huang, S.-L. Tai, B. Singh, P. A. Bogenschutz, X. Zheng, W. Lin, J. Quaas, H. Chepfer, M. A. Brunke, X. Zeng, J. Mulmenstadt, S. Hagos, Z. Zhang, H. Song, X. Liu, M. S. Pritchard, H. Wan, J. Wang, Q. Tang, P. M. Caldwell, J. Fan, L. K. Berg, J. D. Fast, M. A. Taylor, J.-C. Golaz, S. Xie, P. J. Rasch, and L. R. Leung, “Better calibration of cloud parameterizations and subgrid effects increases the fidelity of the E3SM atmosphere model version 1,” Geosci. Model Dev. 15, 2881–2916 (2022).
  * Alizadeh [2022] O. Alizadeh, “Advances and challenges in climate modeling,” Climatic Change 170, 18 (2022).
  * van Nes  _et al._ [2015] E. H. van Nes, M. Scheffer, V. Brovkin, T. M. Lenton, H. Ye, E. Deyle, and G. Sugihara, “Causal feedbacks in climate change,” Nature Climate Change 5, 445–448 (2015).
  * Fischer  _et al._ [1999] H. Fischer, M. Wahlen, J. Smith, D. Mastroianni, and B. Deck, “Ice Core Records of Atmospheric CO2 Around the Last Three Glacial Terminations,” Science 283, 1712–1714 (1999).
  * Stips  _et al._ [2016] A. Stips, D. Macias, C. Coughlan, E. Garcia-Gorriz, and X. S. Liang, “On the causal structure between CO2 and global temperature,” Scientific Reports 6, 21691 (2016).
  * Kaufmann and Juselius [2016] R. K. Kaufmann and K. Juselius, “Testing competing forms of the Milankovitch hypothesis: A multivariate approach,” Paleoceanography 31, 286–297 (2016).
  * Shao and Ditlevsen [2016] Z.-G. Shao and P. D. Ditlevsen, “Contrasting scaling properties of interglacial and glacial climates,” Nat. Commun. 7 (2016), 10.1038/ncomms10951.
  * Ghil and Lucarini [2020] M. Ghil and V. Lucarini, “The physics of climate variability and climate change,” Rev. Mod. Phys. 92, 035002 (2020).
  * Ghil and Simonnet [2020] M. Ghil and E. Simonnet, “Geophysical Fluid Dynamics, Nonautonomous Dynamical Systems, and the Climate Sciences,” in  _Mathematical Approach to Climate Change and its Impacts: MAC2I_, Springer INdAM Series, edited by P. Cannarsa, D. Mansutti, and A. Provenzale (Springer International Publishing, New York, 2020) pp. 3–81.
  * Majda, Timofeyev, and Vanden Eijnden [1999] A. J. Majda, I. Timofeyev, and E. Vanden Eijnden, “Models for stochastic climate prediction,” Proceedings of the National Academy of Sciences of the United States of America 96, 14687–14691 (1999).
  * EPICA community members [2004] EPICA community members, “Eight glacial cycles from an Antarctic ice core,” Nature 429, 623–628 (2004).
  * Zhou and Leung [2010] Y. Zhou and Y. Leung, “Multifractal temporally weighted detrended fluctuation analysis and its application in the analysis of scaling behavior in temperature series,” Journal of Statistical Mechanics: Theory and Experiment 2010, P06021 (2010).
  * Parrenin  _et al._ [2007] F. Parrenin, J.-M. Barnola, J. Beer, T. Blunier, E. Castellano, J. Chappellaz, G. Dreyfus, H. Fischer, S. Fujita, J. Jouzel, K. Kawamura, B. Lemieux-Dudon, L. Loulergue, V. Masson-Delmotte, B. Narcisi, J.-R. Petit, G. Raisbeck, D. Raynaud, U. Ruth, J. Schwander, M. Severi, R. Spahni, J. P. Steffensen, A. Svensson, R. Udisti, C. Waelbroeck, and E. Wolff, “The EDC3 chronology for the EPICA Dome C ice core,” Climate of the Past 3, 485–497 (2007).
  * Akima [1970] H. Akima, “A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures,” Journal of the ACM 17, 589–602 (1970).
  * Miller [2019] J. I. Miller, “Testing Cointegrating Relationships Using Irregular and Non-Contemporaneous Series with an Application to Paleoclimate Data,” Journal of Time Series Analysis 40, 936–950 (2019), _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/jtsa.12469.
  * Kantelhardt  _et al._ [2002] J. W. Kantelhardt, S. A. Zschiegner, E. Koscielny-Bunde, S. Havlin, A. Bunde, and H. E. Stanley, “Multifractal detrended fluctuation analysis of nonstationary time series,” Physica A: Statistical Mechanics and its Applications 316, 87–114 (2002).
  * Agarwal, Del Sordo, and Wettlaufer [2017] S. Agarwal, F. Del Sordo, and J. S. Wettlaufer, “Exoplanetary detection by multifractal spectral analysis,” Astron. J. 153 (2017).
  * Agarwal and Wettlaufer [2022] S. Agarwal and J. S. Wettlaufer, “Minimal data fidelity for stellar feature and companion detection,” Astron. J. 163 (2022).
  * Agarwal, Moon, and Wettlaufer [2012] S. Agarwal, W. Moon, and J. S. Wettlaufer, “Trends, noise and re-entrant long-term persistence in Arctic sea ice,” Proc. Roy. Soc. A 468, 2416–2432 (2012).
  * Moon, Agarwal, and Wettlaufer [2018] W. Moon, S. Agarwal, and J. S. Wettlaufer, “Intrinsic Pink-Noise Multidecadal Global Climate Dynamics Mode,” Phys. Rev. Lett. 121, 108701 (2018).
  * Hasselmann [1976] K. Hasselmann, “Stochastic climate models Part I. Theory,” Tellus 28, 473–485 (1976).
  * Moon and Wettlaufer [2017] W. Moon and J. S. Wettlaufer, “A unified nonlinear stochastic time series analysis for climate science,” Sci. Rep. , 44228 (2017).
  * Moon and Wettlaufer [2019] W. Moon and J. S. Wettlaufer, “Coupling functions in climate,” Phil. Trans. Roy. Soc. A 377, 20190006 (2019).
  * Othmer and Scriven [1971] H. G. Othmer and L. E. Scriven, “Instability and dynamic pattern in cellular networks,” J. Theor. Bio. 32, 507–537 (1971).
  * Levin [1974] S. Levin, “Dispersion and population interactions,” Amer. Natur. 108, 207–228 (1974).
  * Kopell and Howard [1973] N. Kopell and L. N. Howard, “Plane wave solutions to reaction-diffusion equations,” Stud. Appl. Math. 52, 291–328 (1973).
  * Krause  _et al._ [2021] A. L. Krause, E. A. Gaffney, P. K. Maini, and V. Klika, “Modern perspectives on near-equilibrium analysis of Turing systems,” Phil. Trans. R. Soc. A 379, 20200268 (2021).
  * Note [1] We note, however, that a different connotation of stability can be used [37]. Namely, in the context of the longevity of a climate variable, when $b_{12}(t)$ is negative (positive), and thereby provides a weak positive (negative) forcing to variable 1, we can say that variable 2 stabilizes (destabilizes) the presence in the climate system of variable 1.
  * Smale [1976] S. Smale, “A mathematical model of two cells via turing’s equation.” in  _In The Hopf bifurcation and its applications_ , Applied Mathematical Sciences, Vol. 19, edited by J. Marsden and M. McCracken (Springer, New York, NY, 1976) pp. 354–367.
  * Cheng and Redfern [2022] C.-H. Cheng and S. A. T. Redfern, “Impact of interannual and multidecadal trends on methane-climate feedbacks and sensitivity,” Nat. Commun. 13, 3592 (2022).
  * Arneth  _et al._ [2010] A. Arneth, S. P. Harrison, S. Zaehle, K. Tsigaridis, S. Menon, P. J. Bartlein, J. Feichter, A. Korhola, M. Kulmala, D. O’Donnell, G. Schurgers, S. Sorvari, and T. Vesala, “Terrestrial biogeochemical feedbacks in the climate system,” Nat. Geosci. 3, 525–532 (2010).
  * van Groenigen, Osenberg, and Hungate [2011] K. J. van Groenigen, C. W. Osenberg, and B. A. Hungate, “Increased soil emissions of potent greenhouse gases under increased atmospheric CO2,” Nature 475, 214–U121 (2011).
  * Stocker  _et al._ [2013] B. D. Stocker, R. Roth, F. Joos, R. Spahni, M. Steinacher, S. Zaehle, L. Bouwman, Xu-Ri, and I. C. Prentice, “Multiple greenhouse-gas feedbacks from the land biosphere under future climate change scenarios.” Nature Clim. Change 3, 666–672 (2013).
  * Khalil and Rasmussen [1989] M. A. K. Khalil and R. A. Rasmussen, “Climate-induced feedbacks for the global cycles of methane and nitrous oxide,” Tellus 41B, 554–559 (1989).
  * Pierrehumbert [2010] R. T. Pierrehumbert,  _Principles of Planetary Climate_ (Cambridge University Press, Cambridge, 2010).
  * Cuffey and Vimeux [2001] K. Cuffey and F. Vimeux, “Covariation of carbon dioxide and temperature from the vostok ice core after deuterium-excess correction,” Nature 412, 523–527 (2001).
  * Kang and Larsson [2014] J. Kang and R. Larsson, “What is the link between temperature and carbon dioxide levels? a granger causality analysis based on ice core data.” Theor. Appl. Climatol. 116, 537–548 (2014).
  * Persson [2019] E. K. Persson, “What is the link between temperature, carbon dioxide and methane? a multivariate granger causality analysis based on ice core data from dome c in antarctica,” Tech. Rep. U.U.D.M. Project Report 2019:50 (Uppsala University, 2019).
  * Loulergue  _et al._ [2008] L. Loulergue, A. Schilt, R. Spahni, V. Masson-Delmotte, T. Blunier, B. Lemieux, J.-M. Barnola, D. Raynaud, T. F. Stocker, and J. Chappellaz, “Orbital and millennial-scale features of atmospheric ch4 over the past 800,000 years,” Nature 453, 383–386 (2008).