# Probing the timescale dependency of local and global variations in surface
air temperature from climate simulations and reconstructions of the last
millennia

Beatrice Ellerhoff<EMAIL_ADDRESS>Institute of
Environmental Physics, Ruprecht-Karls-Universität Heidelberg, INF 229, 69120
Heidelberg, Germany  Kira Rehfeld  Institute of Environmental Physics,
Ruprecht-Karls-Universität Heidelberg, INF 229, 69120 Heidelberg, Germany  now
at: Geo- und Umweltforschungszentrum (GUZ), Universität Tübingen,
Schnarrenbergstr. 94-96, 72076 Tübingen, Germany

(Published 27 December 2021 in Physical Review E)

###### Abstract

Earth’s climate can be understood as a dynamical system that changes due to
external forcing and internal couplings. Essential climate variables, such as
surface air temperature, describe this dynamics. Our current interglacial, the
Holocene (11,700 yr ago to today), has been characterized by small variations
in global mean temperature prior to anthropogenic warming. However, the
mechanisms and spatiotemporal patterns of fluctuations around this mean,
called temperature variability, are poorly understood despite their socio-
economic relevance for climate change mitigation and adaptation. Here, we
examine discrepancies between temperature variability from model simulations
and paleoclimate reconstructions by categorizing the scaling behavior of local
and global surface air temperature on the timescale of years to centuries. To
this end, we contrast power spectral densities (PSD) and their power-law
scaling using simulated and observation-based temperature series of the last
6000 yr. We further introduce the spectral gain to disentangle the externally
forced and internally generated variability as a function of timescale. It is
based on our estimate of the joint PSD of radiative forcing, which exhibits a
scale break around the period of 7 yr. We find that local temperature series
from paleoclimate reconstructions show a different scaling behavior than
simulated ones, with a tendency towards stronger persistence (i.e.,
correlation between successive values within a time series) on periods of 10
to 200 yr. Conversely, the PSD and spectral gain of global mean temperature
are consistent across data sets. Our results point to the limitation of
climate models to fully represent local temperature statistics over decades to
centuries. By highlighting the key characteristics of temperature variability,
we pave a way to better constrain possible changes in temperature variability
with global warming and assess future climate risks.

Scaling laws of complex systems, fractional Brownian motion, climate
variability, Multiple timescale dynamics, statistical methods, climate
research, linear response theory, power spectral density

††preprint: APS/123-QED

##  I Introduction

The variability of surface air temperature is present on all spatial and
temporal scales, from synoptic and seasonal changes to long-term variations on
periods of years to multi-millennia. On the one hand, it arises from internal
processes, such as the El Niño-Southern Oscillation (ENSO) [1]. On the other
hand, the temperature varies due to external forcing, such as the greenhouse
effect [2, 3]. Understanding the internally generated and externally forced
variability has been suggested to be at least as necessary for evaluating
climate risks for society and ecosystems as projecting the global mean
temperature [4]. Available instrumental observations are limited to a small
time span, leading to challenges in quantifying temperature variability.
Paleoclimate reconstructions extend the characterization of temperature
variability and can be compared to global circulation models (GCMs) [5, 6, 7].
However, discrepancies between model and paleoclimate data remain to be
resolved, especially on the local level and on periods between years and
centuries [8, 9, 10, 11, 12].

Characterizing local temperature variability is crucial for predicting
extremes [6], not only to minimize short-term damage but also to design long-
term strategies, including urban planning and food cultivation [13].
Variability of global temperature on periods above years is relevant to the
understanding of long-term changes [14] as well as climate sensitivity [15].
Assessing the temporal correlation structure of temperature series by means of
scaling behavior and persistence is particularly important for distinguishing
externally forced trends from natural changes [16]. It could affect the
confidence in future projections and attribution studies [17, 18]. Therefore,
one of the main topics to be investigated here is the characteristics of local
and global temperature variability on periods of years to centuries from model
simulations and observation-based data of the last millennia.

To determine how the variability of a temperature series is distributed with
timescales $\tau$, we make use of the power spectral density (PSD) $S(\tau)$,
known as “spectrum”. It can be obtained from the Fourier transform of the
autocorrelation function (see Appendix A) [19, 20]. The spectrum was shown to
often follow a power law

$S(\tau)\sim\tau^{\beta}\,,$ (1)

with spectral exponent $\beta$ and period $\tau$ [21, 22, 23, 24, 25, 26],
especially on decadal-to-centennial scales [27, 28, 29]. We refer to this
behavior (1) as temporal scaling since the temperature signal has no preferred
timescale and is statistically similar across periods $\tau$. The exact
determination of the start and end points of a scaling interval is not part of
this study.

Long-range memory stochastic processes are suitable to describe temperature
signals with temporal scaling [25, 24, 30]. Among those, fractional Gaussian
noise (fGn) is a stationary process and exhibits a spectral exponent
$\beta\in\left(-1,1\right)$ on sufficiently long periods (see Appendix B).
Fractional Brownian motion (fBm) is a nonstationary process that shows
$\beta\in\left(1,3\right)$. The scaling exponent $\beta$ relates to the decay
of the autocovariance function (see Appendix B) and indicates how strongly the
values within a time series are correlated (or anticorrelated). It is
therefore regarded as a measure of the strength of temporal persistence (or
antipersistence) [31, 25].

Figure 1:  Characteristic timescales relevant to surface air temperature (T)
variability of climatic drivers (dark blue) and climate subsystems (yellow)
[32, 33, 34]. The weather and long-term climate is characterized by $\beta>1$
for local and global mean temperature. On interannual to millennial timescales
the statistical properties of temperature fluctuations remain to be
determined, especially at the local scale. The TSI bar highlights the dominant
variations in recent total solar irradiance observations.

Particular scaling behavior with $\beta\approx 2$ [24, 22] is typical for the
weather regime (hours to weeks) and can be explained by atmospheric turbulence
[35, 36]. In the long-term climate, regional and global mean temperatures show
similar behavior ($\beta>1$) [24, 23, 29] due to the presence of nonlinear
processes, such as the temperature-albedo feedback [37]. On timescales between
years to millennia, the temperature is constantly influenced by the
interaction of all climate subsystems and by volcanic, solar, as well as CO2
forcing (Fig. 1). Estimates of the spatially-dependent scaling behavior of
local temperature on these timescales differ [24, 27]. On the global scale,
many studies find $\beta\approx 1$ [28, 29]. However, Lovejoy et al. has
identified a change from the so-called macroweather regime ($\beta\approx 0.8$
on periods of 10 days to 40 yr) to the climate regime ($\beta\approx 1.8$ on
periods from 40 yr to 80 000 yr) [24].

In this manner, previous works find ambiguity in the interpretation of local
and global temperature scaling and it remains to be determined whether
simulations and reconstructions qualitatively agree in scaling behavior
$\beta<1$ or $\beta>1$. The so-called “$1/f-$noise” ($\beta=1$) corresponds to
a process with power spectral density proportional to the period. For
$\beta>1$, the relative contribution

$\frac{\int_{f^{\prime}/2}^{f^{\prime}}S(f)\,\mathrm{d}f}{\int_{f^{\prime}}^{2f^{\prime}}S(f)\,\mathrm{d}f}=\frac{1-2^{\beta-1}}{2^{1-\beta}-1}=2^{\beta-1}$
(2)

to the variance is larger from slow timescales compared to faster ones for all
frequency intervals $f^{\prime}/2$ to $2f^{\prime}$ within a scaling interval
[38]. With increasing $\beta>1$, the fBM is said to exhibit “nonlinear pseudo-
trends” [30] (see Appendix B). Thus, for understanding climate variability and
for modeling purposes, the systematic estimate of the scaling exponent $\beta$
allows to assess the behavior of fluctuation levels across timescales [25].
Moreover, the differentiation between forced and unforced changes poses a
challenge to understanding temperature variability [39, 40]. Beyond the
analysis of Haar fluctuations of a few forcing reconstructions [41, 42, 43],
spectral analysis of climatic drivers and their frequency-dependent linkage to
the temperature response remains incomplete.

We investigate the timescale dependency of local and global surface air
temperature variability by analyzing power spectral densities from a few hours
to a thousand years, thereby extending and improving on earlier work [44, 24,
23, 29]. We use model simulations and observation-based data, which we
introduce in Sec. II. To estimate the PSD and determine its power-law scaling
on periods of 10 to 200 yr, we use state-of-the-art methods described in Sec.
III. This allows us to contrast regional and global spectra (Sec. IV.1),
spatial patterns (Sec. IV.2), and the agreement of simulated and observation-
based estimates (Sec. III.3). Along with that, we discuss the joint PSD from
various radiative forcings, which allows us to calculate the spectral gain and
study the externally forced variability in Sec. III.4. Based on our
reconstruction of the PSD of surface air temperature for the last millennia,
we evaluate the consistency of spectral characteristics across the data sets
considered. In Sec. V, we elaborate on the stronger persistence of temperature
on local than global level as well as remaining discrepancies. Finally, we
discuss how our findings could help improve climate model simulations and
understand Earth’s climate dynamics.

##  II Data

We investigate the timescale-dependent distribution of surface air temperature
variability using model simulations, observation-based data, and radiative
forcing reconstructions. The model simulations include ten transient runs from
GCM experiments [45]. The observation-based data consists of reanalysis data,
instrumental measurements, and the paleoclimate reconstructions from the Past
Global Changes 2k (PAGES2k) network [46]. We use 12 reconstructions of
climatic drivers, including solar, volcanic, orbital, and CO2 forcing. All
temperature and radiative forcing signals are specific to the Mid- and Late-
Holocene (the last 6000 yr), with a focus on the Common Era (0 to 2000 CE).
The supplemental tables S1-S3 [47] summarize their key specifications.

###  II.1 Model simulations

Each of the ten GCM runs considered features a transient, albeit different
forcing and a comparable spatiotemporal resolution. The CESM-LME 1 [48] and
MPI-M LM [49] experiments serve as representative runs of the last millennium.
We analyze them at two temporal resolutions (one month, six hours) to capture
both the high- and low-frequency variability within our available computing
capacities (see Fig. S7 [47]). CESM 1 past 2k [50] is included as a slightly
newer run for the Common Era. To cover the Mid-Holocene, we use simulations
from the IPSL [51] (denoted IPSL-p6k) and ECHAM5/MPI-OM [52] (denoted
ECH5/MPIOM-p6k) of the last 6000 yr. From the TraCE-21k [53] simulation, we
also consider only the last 6000 yr to retain comparability and to avoid
potential biases due to significant shifts in the mean state of climate. The
Mid-Holocene runs were temporally averaged to a bi-monthly resolution to
reduce computational costs. To test for the influence of human-induced climate
change on our results, we include the HadCM3 LM1 simulation [12], covering the
period from 850 to 1850 CE. Furthermore, we compare our results to the pre-
industrial (PI) control runs from CESM-LME 1 and MPI-M LM, as well as the
TraCE-21k-ORB run, which is solely forced by orbital changes.

###  II.2 Observation-based data

In addition to the simulations, we analyze the monthly resolved HadCRUT4
(Hadley Centre/Climatic Research Unit Temperature) instrumental records,
ranging from 1850 to 2019 [54]. However, most of the grid-box time series are
not available as continuous measurements as required for spectral analysis.
Therefore, we retain only those 104 grid boxes with coverage greater than 150
yr after interpolating gaps of up to two months. While the Northern Hemisphere
is comparatively well covered up to 72.5∘N, only nine grid boxes remain for
the Southern Hemisphere. Therefore, this selection comes at the expense of
spatial resolution but offers a higher spectral resolution on longer
timescales. To further explore the potential effect of these spatiotemporal
constraints, we include the ERA5 (European Centre for Medium-Range Weather
Forecasts Reanalysis 5th generation) temperature reanalysis for the years 1979
to 2019 [55]. Along with CESM-LME 1 and MPI-M LM, we analyze the ERA5 data at
both six-hourly and monthly resolution for the same reasons as mentioned
earlier.

In addition to direct temperature observations and reanalysis, we analyze
paleoclimate data. Paleoclimate records hold preserved biological, chemical,
and physical tracers (“proxies”) of past climate. The number of temperature
records from paleoclimate data with sub-centennial resolution is limited.
Recent progress has been made by improved calibration and pseudo-proxy methods
within the PAGES2k network [56]. Therefore, we base our analysis on their
newest global multiproxy database for temperature reconstructions of the
Common Era [46]. It gathers 692 records from trees, ice, sediment, corals,
speleothems, and documentary evidence with a resolution between weeks and
centuries. The records are spread over 648 locations, including all
continental regions and major ocean basins.

For investigating the variability of global mean surface temperature, we use
the seven spatially-weighted statistical reconstructions for the last 2000 yr
provided by PAGES2k [46]. To estimate the mean of local spectra, we choose
records from the PAGES2k database according to their resolution ($\leq 80$yr),
their number of data points ($\geq 20$), their coverage ($\geq 20$ yr), as
well as their maximum hiatuses ($\leq$ 160 yr). To reliably deduce the scaling
of the PSD from individual records, we select the records according to our
scales of interest (Table 1), similar to [57, 25]. Ice core records were
excluded from our analysis since they require additional consideration of
signal-to-noise ratios at the sub-centennial timescales [58, 59].

Table 1:  Requirements on irregularly sampled time series $x(t)$ for analyzing power-law scaling on timescales $\tau\in\left[\tau_{1},\tau_{2}\right]$. We apply this scheme for $\tau_{1}=10$ and $\tau_{2}=200$ yr in Sec. IV.2 and IV.3. Parameter | Value
---|---
Number of data points ($N$) |  $\geq$ 50
Mean temporal resolution ($\langle t_{i+1}$ \- $t_{i}\rangle$) | $\leq\tau_{1}$
Coverage ($t_{N}$ \- $t_{1}$) | $\geq 3\tau_{2}$
Length of hiatuses (max($t_{i+1}$ \- $t_{i}$)) | $\leq 5\tau_{1}$

###  II.3 Radiative forcing

External forcing contributes significantly to temperature variability and is
an essential part of reliable climate projections [40, 60, 61]. We study its
spectral properties using forcing reconstructions, widely implemented in GCM
experiments and coordinated within the Palaeoclimate Model Intercomparison
Project (PMIP3/PMIP4) [62, 63]. This includes five solar [64, 65, 66, 67, 68],
one CO2 [63] and two volcanic [60, 69] forcing reconstructions as well as
Berger’s numerical solution for orbital forcing [70]. Furthermore, we
calculate diurnal insolation changes from the hour angle of the sun [71]. We
also use a more recently published volcanic [72] and high-resolution solar
forcing [73] reconstruction as well as CO2 measurements [74]. We neglect land-
use forcing [75] which is much lower in amplitude and variability than the
other forcings considered here.

All forcing reconstructions are rescaled to radiative forcing equivalents,
which express their respective change in the Earth’s radiation balance in
Watts per square meter ($\mathrm{Wm^{-2}}$). We apply the widely used formula
$5.35\,\mathrm{ln}(\left[\mathrm{CO}_{2}\right]/278\mathrm{ppm})\,\mathrm{Wm^{-2}}$
to rescale CO2 concentrations $\left[\mathrm{CO}_{2}\right]$, given in parts
per million (ppm) [76]. The stratospheric aerosol optical depth (AOD) from
volcanic eruptions is rescaled by $(-20)^{-1}$ $\mathrm{Wm^{-2}}$/AOD [77],
however, the optimal conversion factor is still a matter of debate [78].
Additional uncertainties arise from the wide spread of reconstructions for
volcanic and solar forcing. To account for this and the choice of conversion
factor, we simulate the joint PSD of radiative forcing by a Monte Carlo
approach described in Appendix E. Here, “joint” indicates that the PSD of
radiative forcing is calculated by linear summation of the mean PSD from
different types of climatic drivers, rescaled to their radiative forcing
equivalents.

##  III Methods

Spectral analysis is the primary tool used here for studying the timescale-
dependent variability and scaling of temperature series. To minimize
uncertainties in the spectral analysis of proxy records, we use state-of-the-
art approaches for irregularly sampled time series [79]. Statistical
estimators further test for the agreement between simulations and paleoclimate
data. We apply linear response theory to derive the spectral gain and
investigate the forced temperature response.

###  III.1 Spectral analysis

Power spectral analysis requires the assumption that the underlying time
series can be described as a weakly stationary, stochastic process with time-
independent mean and autocovariance [80]. We therefore linearly detrend all
time series as it is standard for temperature analysis [81, 25, 82, 9]. The
agreement of the PSD from disjoint time intervals in Fig. S13 [47] provides
evidence that stationarity is sufficiently fulfilled. We use the multitaper
method with three windows [83, 84] and chi-square distributed uncertainties to
compute the PSD. The two lowest frequencies were omitted to reduce biases of
the multitaper method [23]. For visual purposes, we apply a logarithmic
Gaussian smoothing filter of constant width (0.005 decibels) [85]. Mean
spectra were calculated by interpolation to the lowest resolution, binning
into equally spaced log-frequency intervals, and taking the average with equal
weights [23]. This requires the statistical independence of the averaged
values [43]. The spectral exponent $\beta$ is calculated by linear regression
to the logarithm of (1) on periods between $\tau_{1}=10$ and $\tau_{2}=200$ yr
after binning the PSD into equally spaced log-frequency intervals to more
uniformly weight the estimate and avoid low-frequency biases [23, 25, 29, 86].
In the case of seven proxy records with an insufficient resolution, the
scaling is estimated on their corresponding spectral resolution, but always at
least between 20 and 200 yr (Fig. S4 [47]). The uncertainty of the spectral
exponent, $\Delta\beta$, is given by the standard error of the linear
regression model $\Delta\beta_{lm}$, except for irregularly temperature
series.

###  III.2  Uncertainties for irregular temperature series

Spectral analysis of proxy records, which are typically not sampled in regular
time steps, is more prone to errors than that of regular time series. We aim
to minimize biases by accounting for the number of data points, temporal
resolution, total coverage, and hiatuses’ length when selecting the records
(Table 1). We find that the mean temporal resolution of a proxy record
approximates well the optimal interpolation time step. Nevertheless, the
interpolation introduces uncertainties which are not captured by
$\Delta\beta_{lm}$. Similar to Laepple et al. [79], we quantify this
additional uncertainty $\Delta\beta_{int}$ in four steps: (1) For each record
with spectral exponent $\beta$, we simulate $N=100$ surrogate time series with
annual resolution and a power-law scaling $\beta_{n}\approx\beta$ and
$n\in[1,N]$. (2) We form the surrogate’s block-average over the proxy record’s
irregular time steps and obtain $N$ surrogate time series at record
resolution. (3) We interpolate the surrogate time series, calculate the
multitaper spectrum, and extract the scaling exponent $\beta_{n,lm}$ from
linear regression in the same way as for the proxy record (Fig. S8 [47]). (4)
We calculate the mean deviation
$\Delta\beta_{int}=\frac{1}{N}\sum_{n=1}^{N}|\beta_{n,lm}-\beta_{n}|$ of the
ensemble. The uncertainty of the individual fits $\Delta\beta_{n,lm}$ is
negligible compared to the mean deviation $\Delta\beta_{int}$. We obtain the
uncertainty of the record’s spectral exponent from both, the uncertainty of
the initial fit $\Delta\beta_{lm}$ and due to interpolation
$\Delta\beta_{int}$ via quadratic summation:
$\Delta\beta=\sqrt{(\Delta\beta_{lm})^{2}+(\Delta\beta_{int})^{2}}$.

###  III.3 Statistical analysis of spectral exponents

We quantify the agreement of simulated and reconstructed $\beta$-values using
percent agreement, categorical agreement, and Kappa statistics. Beforehand, we
extract the simulated temperature at the proxy record location by bilinear
interpolation of neighboring grid boxes to achieve the best possible
comparability between record and simulation. Percent agreement $p_{0}$ gives
the percentage of locations at which the confidence range
$\beta\pm\Delta\beta$ from simulation and reconstruction overlap. The
agreement by category, here referred to as categorical agreement $p_{c}$, is
calculated with the help of $\nu=0.32$, the mean uncertainty of $\beta$ from
all proxy records considered. We then assign the three categories low
($\beta<1-\nu$), high ($1+\nu\leq\beta$) and intermediate
($1-\nu\leq\beta<1+\nu$) to the spectral exponent $\beta$. The intermediate
regime prevents incorrect assignment. To verify the reliability of categorical
agreement, we calculate the Kappa statistics

$\kappa=(p_{c}-p_{e})/(1-p_{e})$ (3)

with expected percent agreement $p_{e}$ by category [87]. The latter can be
obtained from $p_{e}=\frac{1}{N^{2}}\sum^{3}_{c=1}n_{c,m}n_{c,p}$ where $c$ is
the category, $N$ the number of locations and $n$ the number of times that
models ($m$) and proxy records ($p$) have predicted category $c$. The
$\kappa$-coefficient quantifies the reliability from no agreement beyond
chance ($\kappa=0$) to full agreement ($\kappa=1$). Negative $\kappa$
indicates agreement that is beyond change, for example, due to systematic
biases.

###  III.4 Spectral gain

We investigate how climatic drivers influence the global mean temperature at
period $\tau$ by calculating the spectral gain

$G^{2}(\tau)=\dfrac{S_{T}(\tau)}{S_{F}(\tau)}\,.$ (4)

Here, $S_{T}(\tau)$ is the PSD of the global mean temperature and
$S_{F}(\tau)$ the PSD of radiative forcing (see also Appendix C). The gain
requires the assumption that the global mean temperature can be well
approximated as a linear function of the forcing [88, 89, 26] and that
different types of radiative forcing add linearly [90, 91, 92, 93]. To this
end, we focus on timescales between years and centuries when additivity is a
valid assumption and nonlinearities in the global mean temperature are
sufficiently small [42, 43]. The main practical problem that confronts us is
that the gain might be subject to a sampling bias due to our data sets choice.
Therefore, we perform a Monte Carlo simulation of the PSD of radiative forcing
and the global mean temperature, as well as the spectral gain as described in
Appendix E.

##  IV  Results and discussion

###  IV.1  Global mean and mean of local spectra

Figure 2:  (a) Mean power spectral densities (PSD) of local temperature from
model simulations and observation-based data on periods from hours to 1000 yr
for the Holocene. (b) PSD of global mean temperature. The dashed lines with
slope $\beta$ and arbitrary $y$-intercept in the log-log graph indicate the
scaling behavior for visual comparison. The ensemble means (black solid lines)
were formed using equal weights across the model group M0 (see Table S1 [47]).

In order to study the timescale dependency of global mean temperature, we
present its power spectral density in Fig. 2 (b). It shows the characteristic
background continuum, spectral peaks, and higher harmonics associated with the
diurnal and annual cycle. Overall, the PSDs tend to agree between the data
sets, albeit with some differences on the interannual scale and when compared
to the Trace21k ORB run. The Trace21k-ORB run is solely forced by orbital
changes and therefore shows less variability than the ensemble mean. The broad
spectral peak on interannual periods reveals an artificially amplified ENSO in
the shared MPI-M LM and ECHAM5/MPI-OM ocean component [94]. For a better
visibility, PI control runs are separately shown in the supplementary Fig. S6
[47]. Overall, the PSD largely agrees among different data sets, especially
towards shorter timescales.

We find a power-law scaling of $\beta\approx 1$ on timescales longer than 10
yr in line with previous results [24, 28, 25]. The PSD decreases more strongly
towards shorter periods, which is characteristic of the weather regime [24,
36]. Similar to Nilsen et al. [25], we find no evidence for significant
changes in scaling behavior around the centennial scale. One limitation of
previous work that found scale breaks is that the spectra were estimated
across nonstationary shifts in climate, such as the deglaciation [29], and
with a change in proxies and archives [23].

We present the area-weighted mean spectra of the local (grid box) temperature
in Fig. 2(a). Compared to the global mean in Fig. 2(b), the power increases
and the spectral slope decreases, in line with [81]. The spectra agree on
periods below 10 yr, except for the artificially amplified ENSO signal
mentioned earlier. Moreover, we find a narrow peak at 13 yr, associated with
an unrealistic variability in the northern North Atlantic of the TraCE-21k
run, similar to [95, 96]. Remarkably, the decadal-to-centennial variability of
the reconstructed temperature is increased by one to two orders of magnitude
compared to the simulations. The spectral exponent is smaller for models
($\beta<1$) compared to paleoclimate data ($\beta\approx 1$).

This finding verifies that models show less regional temperature variability
and that the mismatch increases towards longer timescales. The results are
robust to sampling from the PAGES2k database and the influence of
anthropogenic climate change (Fig. S10 [47]). One shortcoming of forming the
area-weighted mean PSD is that the uncertainty quantification requires the
assumption of independent spatial degrees of freedom of the temperature field.
Due to the presence of spatial correlations, an estimate of the effective
spatial degrees of freedom and their dependence on the underlying timescale
would be needed to resolve this limitation [97].

###  IV.2 Spatial patterns of persistence

Figure 3:  Local temperature persistence on timescales from $\tau_{1}=10$ to
$\tau_{2}=200$ yr across multiple climate simulations and selected proxy
records from the PAGES2k database. Colors from blue to red indicate the
scaling behavior ranging from $\beta=-1$ to $\beta=3$. Symbols indicate the
scaling of proxy records from different natural archives. The background of
each panel shows the $\beta$-values fitted to the PSD of the local grid box
temperature from simulations. Zonal mean values (dashed curves) are given next
to the map, with means (solid curves) over latitude-intervals (with breaks at
$-60$, $-30$, $0$, $30$, and $60^{\circ}$ N) and gray shaded confidence
intervals. The spatial coverage of proxy records is not sufficient for robust
mean estimates, which is why only simulation data are shown here.

To further investigate the mismatch on local scaling properties, we compare
the spatial dependence of temperature persistence from simulations and
paleoclimate data in Fig. 3. The simulations largely exhibit small-magnitude
scaling exponents ($-1<\beta<1$), whereas proxy records were found to also
show $\beta>1$. In this manner, the magnitude of local temperature
fluctuations from model simulations often shows no dependence on the decadal-
to-centennial timescale. However, approximately half of the proxy records show
a variance that grows on increasingly long periods (see also Fig. S11 [47]).

From both simulations and paleoclimate data, we can strengthen the argument by
Fredriksen et al. [81] that there is no latitudinal dependence of $\beta$
(Fig. 3), in contrast to previous studies, suggesting a possible linkage to
the strength of the seasonal cycle [23]. Inspecting the simulations’
$\beta$-values (background of Fig. 3), we find a small land-sea contrast.
Strongest scaling occurs in the Southern Oceans in line with previous findings
[81]. Ocean-sea ice interactions with characteristic timescales of the order
of centuries and a generally increased internal variability over the oceans
might explain these results.

We find generally lower values for the slope $\beta$ in the ENSO and Indo-
Pacific region. This could be attributed to the fact that (quasi-)oscillatory
signals, such as active modes of internal variability, are reflected in the
PSD as broad peaks and hence cannot be described by a scaling law. On the
other hand, this finding is stronger in PI control runs compared to fully
forced runs [Fig. 3 (c-f)]. Thus, residual effects of the recent global
warming trend might play an additional role [98]. A systematic bias becomes
clear from the spatially almost uniform $\beta$-values of Trace21k-ORB [Fig. 3
(h)]. In line with Fig. 2, we explain this by the lack of forcing mechanisms
on interannual to multi-decadal timescales in the aforementioned simulation.

Marine and lake sediments, as well as the archived documents, follow the
general trend of increased $\beta$-values compared to simulations. Tree ring
records agree well with most simulations in North America and Siberia, but not
necessarily at the coast of Australia and northern Europe. Discrepancies such
as those in southern South America could reflect the proxies’ strength in
representing local conditions, for example, topography. However, noise sources
in the climate signal recording and preservation, such as bioturbation, can
influence proxy records. Further separating the signal content from noise
sources in paleoclimate reconstructions can help refine our findings [99,
100].

###  IV.3 Statistical agreement of temperature persistence

Figure 4:  Percentage agreement $p_{0}$, categorical agreement $p_{c}$ and
inter-rater reliability $\kappa$ of local temperature persistence from
simulations and paleoclimate data. The measures were calculated from a set of
bilinearly interpolated simulation records and the proxy record at 23
different locations. Missing orange bars indicate no agreement beyond chance
and, therefore, zero inter-rater reliability ($\kappa=0$).

We further investigate the question of temperature scaling by a statistical
analysis of $\beta$-values from simulations and paleoclimate reconstructions.
It is based on the detailed uncertainty quantification outlined in subsection
III.2. Our results show that reconstructions and simulations agree in less
than 30% of locations within the scope of uncertainties (Fig. 4). To single
out the scaling behavior of temperature signals, we study the agreement by
category. We find approximately 25% of agreement within the categories
$\beta<1-\nu$ (low) and $\beta>1+\nu$ (high). Although widely accepted [101],
categorical and percentage agreement suffer from the limitation to ignore any
agreement by chance. Therefore, we investigate the $\kappa$-statistics (orange
bar in Fig. 4) and verify that there is no agreement beyond chance
($\kappa=0$) for almost all models. Only MPI-M LM and HadCM3 LM1 show any, if
poor agreement ($\kappa\approx 0.1$), whereas Trace21k-ORB shows even lower
agreement than expected by chance ($\kappa<0$) due to its systematic bias.

The disagreement could be attributed to both paleoclimate data and
simulations. A systematic bias could arise, for example, through the recent,
nonstationary global warming trend. Therefore, we repeat our analysis with all
time series cut at 1850. In particular, anthropogenic warming slightly
increases long-term temperature variability and thus scaling behavior, but not
significantly (Figures S6, S9, and S10 [47]). Further uncertainties could
arise from our choice of statistical estimator for the scaling exponent
$\beta$. Maximum likelihood estimation (MLE) should generally be preferred
over linear regression (LR) because of its mathematical soundness and
skillfulness [102]. We find that MLE is indeed more accurate for regular time
series with $\beta>0$ (Fig. S14 [47]). However, LR allows for estimation of
$\beta<1$, unlike MLE which assumes $\beta>1$ [102]. In addition, for the
characteristics of our empirical data, the differences between the two methods
are not significant for $\beta>0$ (Fig. S15 [47]). Therefore, linear
regression represents the preferred estimator for our analysis. Regardless of
the chosen estimator, we observe a slight tendency towards increased scaling
exponents for irregularly sampled data (Fig. S15 [47]), similar to Lucke et
al. [100]. Our uncertainty quantification carefully accounts for these
potential errors due to irregular sampling and interpolation by simulating
their influence using surrogates (Fig. S8 [47]).

We do not expect other systematic biases for the paleoclimate data since we
base our results on multiple archives and proxies, and no systematic spatial
pattern is discernible (Fig. S11 [47]). In particular, the cross-correlations
between the 23 proxy datasets are weakly positive (0.02 on average with 95%
quantiles of -0.17 to 0.21). The assumption of spatial independence necessary
for robust statistical analysis (Fig. S16 [47]) therefore appears fully
satisfied. The models’ resolutions are another possible element of uncertainty
that impacts variability over a wide range of timescales [103, 104, 105]. We
here facilitate inter-model comparison by using state-of-the-art GCMs with
comparable spatial and temporal resolutions, but computational costs precluded
higher resolutions. The latter might be necessary to improve the
representation of decadal variability and response to external forcing. In
particular, the increased scaling exponents ($\beta>1$) from paleoclimate data
could indicate that nonlinear processes from an interactive carbon cycle and
dynamical ice sheets might not be sufficiently represented in models.

###  IV.4 The forced temperature response

Figure 5:  Power spectral densities from radiative forcings. Details on the
reconstructions considered here are summarized in Table S2 and Figure S5 [47].

Climatic drivers are not constant in time and thus affect the surface air
temperature on multiple timescales. To investigate the forced temperature
response, we present spectra for the main climatic drivers in Fig. 5. The PSD
of orbital forcing consists of the diurnal and annual cycle as well as a
background continuum on longer timescales. Higher harmonics on monthly
timescales were omitted. We calculate the mean volcanic, solar, and CO2
spectra using an equally weighted average of spectra from multiple data sets
(Fig. S5 [47]). The CO2 forcing follows the orbital forcing. The PSD of solar
forcing again contains more power and has a pronounced peak around the 11 yr
solar cycle. Multiple theories and paleoclimate reconstructions suggest the
increased variability on centennial to millennial periods due to the long-term
behavior of solar activity [106].

Volcanic forcing dominates interannual to centennial scales and undergoes a
scale break around the period of 7 yr, estimated using the goodness of fit
[102]. Above decadal scales, it follows a white noise spectrum with constant
variance. However, the intermittency of volcanic eruptions might have led to
biases in the spectral characteristics [42]. We verify our results using an
analytical approach described in Appendix D. Remarkably, the derived PSD of an
ideal, intermittent time series with Poisson distributed return times explains
our findings. We further demonstrate the scale break by a Monte Carlo
simulation of the joint PSD of radiative forcing in Fig. 6 (a). This finding
raises the question of how the spectrum with a scale break translates into the
continuous spectrum in Fig. 2 (b).

We address this question by calculating the spectral gain (4) on periods
between years and centuries in Fig. 6 (b). Here, observation-based data
include HadCRUT4, ERA5, and PAGES2k again. To account for the model artifacts
explained above, we calculate the gain from the model simulation group M0 and
together with group M+ (see Table S1 [47]). We find that the spectral gain is
similar from observation-based data and the model simulation group M0, which
is the one without artificially amplified ENSO. This suggests that both follow
a similar distribution of timescale-dependent variability, as already
indicated by Fig. 2(b). Large parts of the gain show constant behavior, which
is most pronounced in M0. In a simplified way, the gain might be approximated
by an ideal linear amplifier or damper of the forcing with comparable internal
variability on all timescales. However, we also find a dip around decadal
scales, which is strongest in the gain from measurements. Inspecting Fig. 6
(a), this can be explained by forming the ratio between a spectrum with a
scale break ($\beta>1\rightarrow\beta\approx 0$) and one with moderate scaling
($\beta\approx 1$).

From this standpoint, internal variability slightly grows on periods from
years to centuries when slow processes in the oceans, vegetation, land
surface, and cryosphere become increasingly active (Fig. 1). While the model
simulations follow this general pattern, they may not represent its amplitude
correctly, for example, due to the lack of feedback mechanisms. In addition, a
too high model diffusivity could cause the suppression of low-frequency
variability in model simulations due to a faster energy dissipation over
temporal scales [8]. The PAGES2k multiproxy reconstruction, stemming from
palaeoclimate data, possibly underestimates internal variability on
interannual scales. However, the mean variance ratios in Fig. 5 (b) of the
model estimates agree with those from observations in the global mean. This
leaves us with a conundrum: the global mean temperature based on model
simulations and observations is mostly consistent in its variability, scaling,
and response to forcing. Notwithstanding, locally, the models show a much
lower variance on longer timescales and different scaling behavior than
reconstructions. Thus, it appears that the statistics of local fluctuations
need to be optimized in models but without significantly altering global
properties. To this end, the study of unforced (“spontaneous”) oscillations
[107] and abrupt transitions [108, 109] in the climate system is one promising
approach to improve the representation of local variation. Furthermore,
higher-resolved ocean and atmosphere models with additional mechanisms such as
ice sheet dynamics and an interactive carbon cycle might increase long-range
dependence and persistence of local temperature in the future.

Figure 6:  Monte Carlo simulation of PSD (a) and spectral gain (b) using
temperature and forcing reconstructions as well as model simulations. Shaded
confidence intervals lie between the 5% and 95 % quantiles. We consider only
models from the groups M0 and M+ (Table S1 [47]) to exclude model artifacts
and to represent the historical temperature response in the best possible way.
Notably, M+ contains those simulations with amplified ENSO [94]. (b) Dashed
lines indicate the mean variance ratio
$\langle{S_{T}}\rangle/\langle{S_{F}}\rangle$.

##  V Conclusion

In summary, we have investigated the question of temperature variability on
the timescale of years to centuries. To this end, we have presented power
spectral densities for both local and global surface air temperature from
simulation and observation-based data of the last millennia. On this basis, we
concluded that locally there is a stronger scaling and increased variance in
reconstructions as compared to simulations. Using statistical analysis, we
found that local temperature series extracted from simulations and
paleoclimate reconstructions show different scaling behavior, with proxy
records hinting at a stronger persistence. Furthermore, we have largely
extended the spectral analysis of climatic drivers by estimating the joint PSD
from CO2, solar, volcanic, and orbital forcing using Monte Carlo simulation.
Hereby, we discovered a scale break at the period of approximately 7 yr.
Moreover, we have presented the spectral gain, describing the timescale-
dependent forced temperature response. We found that it is mostly consistent
across data sets and indicates an increasing internal variability on
timescales of decades to centuries.

Our analysis of the spectral gain was limited to global average values and
those timescales where linearity can be reasonably assumed [43, 42, 110].
Nonlinearities are inherent to the climate system, for example, due to the
temperature-albedo feedback. Thus, it will be necessary to examine their
possible effects on multiple spatiotemporal scales to further extend this
work. Studying nonlinearities could also shine new light on the mechanisms of
scaling in Earth’s climate, which are not yet fully understood and might be
linked to nonlinearities as well [6]. Furthermore, we have focused on the
current interglacial, the Holocene. This is because climate variability has
been demonstrated to depend on the mean climate state [82]. Furthermore, major
shifts in climate could potentially violate the basic assumption of weak
stationarity for spectral analysis. Thus, the conclusions laid out here cannot
be readily applied to other climate states, such as glacial periods, which is
an issue for future studies. Clearly, understanding the dependence of
temperature variability on global warming demands additional work.

Ideally, our findings should be replicated by employing models with increased
internal variability on longer timescales and paleoclimate data that provides
improved spatiotemporal resolution. In particular, investigating the
relationship between spatial and temporal disagreement is a key task for
future analyses. Optimized analysis of noise sources and spectral analysis of
(pseudo-)proxy records could help to expand the data basis of proxy records
with decadal resolution [111, 112, 59]. Regarding climate models, an improved
representation of processes that increase Earth’s long-term memory, such as an
interactive carbon cycle and dynamical ice sheets, might strengthen the long-
range dependence and persistence of surface air temperature. A better
understanding of unforced low-frequency oscillations as well as abrupt changes
will be necessary to improve the representation of local fluctuations and
could further help to understand nonlinear feedback and possible bifurcations
in the climate system. Future studies could also continue to explore how
internally generated and externally forced variability compares on different
spatial scales. Research on the interrelation between internal and forced
changes, as well as local, regional, and global variability, might prove
important and could be conducted using single-forcing experiments from
ensembles of model simulations.

Managing climate risks requires a detailed understanding of temperature
variability. Locally and on timescales between years and centuries, there is
an urgency to address discrepancies to make further progress in climate
modeling. In this study, we have singled out the key characteristics of
temperature variability and showed that the timescale dependency of local
temperature variations from observation-based data and model simulations
differs. Our results have demonstrated that the scaling behavior and spectral
gain are easy-to-use yet effective and promising tools for investigating
variability in Earth’s dynamic climate.

Code to reproduce all figures is available at [113].

###### Acknowledgements.

This manuscript is based upon data provided by the World Climate Research
Programme’s Working Group on Coupled Modelling, which is responsible for CMIP
and PMIP. We thank the research groups listed in Tables S1 and S2 for
producing and making available their data from model outputs, measurements,
paleoclimate, and forcing reconstructions. This study benefited from
discussions within the CVAS working group, a working group of the Past Global
Changes (PAGES) project. We thank T. Gasenzer, T. Kunz, and N. Weitzel for
discussions and J. Bühler, M. Casado, M. Schillinger, and E. Ziegler for
helpful comments on the manuscript. We are grateful to Aimé Fournier and one
anonymous referee for their constructive and valuable review. This research
has been funded by the Heidelberg Graduate School for Physics, by the PalMod
project (subproject no. 01LP1926C), and by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) – project no. 395588486.

##  Appendix A  Relation between power spectral density and variance

The power spectral density of a weakly-stationary, stochastic process is given
by the Fourier transform of the autocorrelation $S(f)=\mathcal{F}\\{R(h)\\}$
with frequency $f$ and lag $h=t_{2}-t_{1}$ between two points in time [19,
20]. For zero lag and zero mean, the integral of the PSD corresponds to the
variance of the signal [80]. Instead of frequency, we use the period
$\tau=1/f$ to express the PSD and spectral gain. The integration of expression
(1) is divergent for $\beta<1$ and $f\rightarrow\infty$ which requires a high-
frequency cut-off, such as described by Lovejoy et al. [38]. In case of
temperature time series considered here, this is naturally defined by the
temporal resolution, setting the maximum frequency.

##  Appendix B  Autocovariance of long-range memory processes

Fractional Brownian motion (fBm) and fractional Gaussian noise (fGn) are fully
described by their correlation properties [30, 114], summarized below. The
autocovariance function of fBm $B(t)$ reads

$\displaystyle\gamma(t^{\prime},t)=\langle B(t^{\prime})B(t)\rangle$
$\displaystyle=\frac{V_{\beta}}{2}\left(|t|^{\beta-1}+|t^{\prime}|^{\beta-1}-|t^{\prime}-t|^{\beta-1}\right)$
$\displaystyle\propto
1+\Big{|}\frac{t^{\prime}}{t}\Big{|}^{\beta-1}-\Big{|}1-\frac{t^{\prime}}{t}\Big{|}^{\beta-1}$
(5)

for $1<\beta<3$. $V_{\beta}$ is a positive constant factor related to
$\langle(B(t^{\prime})-B(t))^{2}\rangle=V_{\beta}\,|t^{\prime}-t|^{\beta-1}$.
By definition, fGn is the series of stationary increments $B(t^{\prime})-B(t)$
and shows spectral exponent $-1<\beta^{\prime}=\beta-2<1$ for $f\ll
1/\pi\Delta t$ with $\Delta t=t^{\prime}-t$. Its autocovariance

$\displaystyle\gamma(h)$
$\displaystyle=\langle(B(t+1+h)-B(t+h))(B(t+1)-B(t))\rangle$
$\displaystyle=\frac{V_{\beta}}{2}|h-1|^{\beta^{\prime}+1}-2|h|^{\beta^{\prime}+1}+|h+1|^{\beta^{\prime}+1}\,,$
(6)

depends only on the lag $h\in\mathbb{Z}$, where we set $\Delta t=1$ without
loss of generality. The fGn has a power spectrum of the form [115]

$S(f)\propto\frac{\sin^{2}(\pi\Delta tf)}{|2\pi\Delta
tf|^{\beta^{\prime}+2}}\,,$ (7)

with the slowly varying factor

$\sin^{2}(\pi\Delta
tf)\underset{f/f_{\mathrm{max}}\rightarrow\,0}{\xrightarrow{\hskip
34.14322pt}}(\pi\Delta t)^{2}f^{2},\qquad f_{\mathrm{max}}=1/\pi\Delta t\,.$
(8)

Considering positive frequencies $f>0$, the spectrum (7) can be approximated
by the power law $S(f)\sim 1/f^{\beta^{\prime}}$ if $f\ll f_{\mathrm{max}}$.
For $f\gtrsim f_{\mathrm{max}}$, however, the fGn has a similar spectral shape
to fBm [114]. We account for this by considering sufficiently long periods. To
give an example, $10^{0.58}$ yr${}^{-1}\lesssim f_{\mathrm{max}}\lesssim
10^{2.7}$ yr-1 corresponds to $6$ h $\lesssim\,\Delta t\,\lesssim 1$ mon.

For all $|t^{\prime}/t|\gg 1$, the covariances (B) keep growing for $\beta>2$
(persistence) and stay bounded for $\beta<2$ (antipersistence). As a result,
equation (B) involves “nonlinear pseudo-trends” [30] for $B(t^{\prime})$
conditioned on $B(t)$, which diverge for $\beta>2$ and converge for $\beta<2$.
According to Eq. (B), fGn is persistent for $\beta^{\prime}>0$ and
antipersistent for $\beta^{\prime}<0$. Ordinary Brownian motion corresponds to
$\beta=2$ and white noise to $\beta^{\prime}=0$. The sequence of partial sums
of the autocovariance function diverges for fGn with $\beta^{\prime}>0$ and
fBm with $\beta>2$. The process is nonsummable and said to possess long-range
memory.

##  Appendix C  Spectral gain for linear systems

In a time-invariant linear system, the output

$y(t)=\int_{-\infty}^{\infty}h(u)x(t-u)\mathrm{d}u$ (9)

is given by the input time series $x(t)$ and the impulse response function
$h(u)$ [80]. The Fourier transform
$H(f)=\mathcal{F}\\{h(u)\\}=G(f)e^{i\phi(f)}$ gives the frequency response
function, also called the transfer function. $G(f)$ and $\phi(f)$ are the gain
and phase, respectively. The integral (9) corresponds to a product in
frequency space $\mathcal{F}\\{y(t)\\}=H(f)\mathcal{F}\\{x(t)\\}$. This
relates the PSD of the output $S_{y}(f)$ to the one of the input $S_{x}(f)$
via

$S_{y}(f)=|H(f)|^{2}S_{x}(f)=G^{2}(f)S_{x}(f)\,.$ (10)

##  Appendix D  Analytical solution to the PSD of intermittent volcanic
forcing

We investigate the power spectral density of intermittent volcanic forcing by
approximating the eruption time series in a simplified way as a stochastic
signal $X(t)=\delta(t-t_{i})$. This function is zero at all times except
$t_{i}$, when an event of unique amplitude occurs. We denote
$T_{i}=t_{i}-t_{i-1}$ the time intervals between two events. We use the fact
that the PSD cannot be calculated only from the covariance, but also from the
Laplace transform $S(X,f)=2\lim_{\epsilon\to
0}\langle|\mathcal{L}(X(t),\frac{\epsilon}{2}-2\pi if)|^{2}\rangle$[116].
Based on this approach, the power spectral density

$S(f)=\mu_{T}\dfrac{1-|\rho(f)|^{2}}{|1-\rho(f)|^{2}}\,,\,f>0$ (11)

becomes a function of the Fourier transform of the probability density
function $\rho(f)=\mathcal{F}\\{\rho(T)\\}$ and the inverse mean interval
between two events $\mu_{T}=\langle T\rangle^{-1}$ [116, 117]. An
exponentially decaying probability distribution
$\rho(T)=\mu_{T}\,\mathrm{exp}(-\mu_{T}T)\Theta(T)$ for volcanic forcing is
suggested [118], and we have checked this for the data sets considered. The
Fourier transform reads $\rho(f)=\mu_{T}(\mu_{T}+2\pi if)^{-1}$ such that
$1-|\rho(f)|^{2}=|1-\rho(f)|^{2}$. As a consequence, the PSD (11) takes a
constant value. We can observe this white noise behavior in Fig. 5 and 6 (a)
on timescales longer than a few years, which is on the order of characteristic
return times for eruptions. Below these timescales, the variability
considerably drops. This analytical result provides an independent
verification of the PSD for volcanic forcing and its scale break.

##  Appendix E  Monte Carlo sampling of the spectral gain

We simulate the spectral gain (4), as well as the PSD of global mean
temperature and the joint PSD of radiative forcing using a Monte Carlo
approach with $N=1000$ realizations to account for sampling biases. The PSD of
global mean temperature is sampled for three groups: the observation-based
data, the model simulations from group M0, and those from M0 together with M+
(Table S1 [47]). Here, only models from the groups M0 and M+ are considered to
exclude model artifacts and to represent the historical temperature response
in the best possible way.

We sample the simulation-based PSD from the average PSD of the simulations
using uniformly distributed random weights. To obtain the observation-based
PSD, we use the global mean temperature from HadCRUT4, ERA5, and a 7000-member
reconstruction ensemble provided by PAGES2k [46]. This ensemble allows us to
sample the PSD by randomly selecting one ensemble member and form the mean of
its spectrum with that of the ERA5 and HadCRUT4 temperature. The joint PSD of
radiative forcing is calculated from all forcing reconstructions considered in
this work except the Fröhlich et al. solar forcing, which has too low temporal
resolution above interannual scales (Table S3 and Fig. S5 [47]). We assume the
PSD of CO2 and orbital forcing as fixed since its spectral power is
comparatively low on multi-decadal scales. We sample the PSD of solar forcing
by using uniformly distributed weights when forming the average PSD of all
solar reconstructions. Similarly, the PSD of volcanic forcing is obtained. In
addition, we randomly vary the conversion factor between $(-18)^{-1}$ and
$(-25)^{-1}$ $\mathrm{Wm^{-2}}$/AOD [78]. The joint PSD of radiative forcing
is calculated by linear summation of the PSD from CO2, orbital, solar, and
volcanic forcing.

Using this sampling scheme, our Monte Carlo produces two outcomes: First, we
compute the PSD of global mean temperature and the joint PSD of radiative
forcing by simulating an ensemble of $N$ realizations for both forcing and
response. Second, we sample the spectral gain directly from the quotient (4)
in each of the $N$ realizations. In both cases, the average of the generated
$N$-member ensemble and its 5% and 95% quantiles constitute the result of our
Monte Carlo simulation.

## References

  * Bjerknes [1966] J. Bjerknes, A possible response of the atmospheric Hadley circulation to equatorial anomalies of ocean temperature, Tellus 18, 820 (1966).
  * Arrhenius [1896] S. Arrhenius, XXXI. On the influence of carbonic acid in the air upon the temperature of the ground, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 41, 237 (1896).
  * Fourier [1824] J. B. J. Fourier, Remarques générales sur les températures du globe terrestre et des espaces planétaires, Annales de Chimie et de Physique , 136 (1824).
  * Katz and Brown [1992] R. W. Katz and B. G. Brown, Extreme events in a changing climate: Variability is more important than averages, Climatic Change 21, 289 (1992).
  * Ghil and Lucarini [2020] M. Ghil and V. Lucarini, The physics of climate variability and climate change, Reviews of Modern Physics 92, 10.1103/RevModPhys.92.035002 (2020).
  * Franzke et al. [2020] C. L. E. Franzke et al., The Structure of Climate Variability Across Scales, Reviews of Geophysics 10.1029/2019RG000657 (2020).
  * Tierney et al. [2020] J. E. Tierney et al., Past climates inform our future, Science 370, 10.1126/science.aay3701 (2020).
  * Laepple and Huybers [2014a] T. Laepple and P. Huybers, Global and regional variability in marine surface temperatures, Geophysical Research Letters 10.1002/2014GL059345 (2014a).
  * Laepple and Huybers [2014b] T. Laepple and P. Huybers, Ocean surface temperature variability: Large model-data differences at decadal and longer periods, Proceedings of the National Academy of Sciences of the United States of America 111, 16682 (2014b).
  * Parsons  _et al._ [2017] L. A. Parsons, G. R. Loope, J. T. Overpeck, T. R. Ault, R. Stouffer, and J. E. Cole, Temperature and precipitation variance in CMIP5 simulations and paleoclimate records of the last millennium, Journal of Climate 10.1175/JCLI-D-16-0863.1 (2017).
  * Ljungqvist  _et al._ [2019] F. C. Ljungqvist, Q. Zhang, G. Brattström, P. J. Krusic, A. Seim, Q. Li, Q. Zhang, and A. Moberg, Centennial-Scale Temperature Change in Last Millennium Simulations and Proxy-Based Reconstructions, Journal of Climate 32, 2441 (2019).
  * Bühler  _et al._ [2021] J. C. Bühler, C. Roesch, M. Kirschner, L. Sime, M. D. Holloway, and K. Rehfeld, Comparison of the oxygen isotope signatures in speleothem records and ihadcm3 model simulations for the last millennium, Climate of the Past 17, 985 (2021).
  * Anderson  _et al._ [2019] W. B. Anderson, R. Seager, W. Baethgen, M. Cane, and L. You, Synchronous crop failures and climate-forced production variability, Science Advances 5, 10.1126/sciadv.aaw1976 (2019).
  * Crucifix  _et al._ [2017] M. Crucifix, A. D. Vernal, and C. Franzke, Centennial To Millennial Climate Variability, Past Global Changes Mag. 25, doi:10.22498/pages.25.3 (2017).
  * Rypdal  _et al._ [2018] M. Rypdal, H. B. Fredriksen, E. Myrvoll-Nilsen, K. Rypdal, and S. H. Sørbye, Emergent scale invariance and climate sensitivity, Climate 6, 1 (2018).
  * Franzke and O’Kane [2017] C. L. E. Franzke and T. J. O’Kane, eds.,  _Nonlinear and Stochastic Climate Dynamics_ (Cambridge University Press, Cambridge, 2017).
  * Barnett et al. [1999] T. P. Barnett et al., Detection and attribution of recent climate change: A status report, Bulletin of the American Meteorological Society 80, 10.1175/1520-0477(1999)080¡2631:DAAORC¿2.0.CO;2 (1999).
  * Bindoff et al. [2013] N. Bindoff et al., Detection and attribution of climate change: From global to regional, in  _Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change_, edited by I. P. o. C. Change (Cambridge University Press, 2013).
  * Wiener [1930] N. Wiener, Generalized harmonic analysis, Acta Mathematica 55, 117 (1930).
  * Khintchine [1934] A. Khintchine, Korrelationstheorie der stationären stochastischen Prozesse, Mathematische Annalen 109, 604 (1934).
  * Wunsch [2003] C. Wunsch, The spectral description of climate change including the 100 ky energy, Climate Dynamics 20, 353 (2003).
  * Franzke [2012] C. Franzke, Nonlinear trends, long-range dependence, and climate noise properties of surface temperature, Journal of Climate 10.1175/JCLI-D-11-00293.1 (2012).
  * Huybers and Curry [2006] P. Huybers and W. Curry, Links between annual, Milankovitch and continuum temperature variability, Nature 441, 329 (2006).
  * Lovejoy [2015] S. Lovejoy, A voyage through scales, a missing quadrillion and why the climate is not what you expect, Climate Dynamics 44, 3187 (2015).
  * Nilsen  _et al._ [2016] T. Nilsen, K. Rypdal, and H. B. Fredriksen, Are there multiple scaling regimes in Holocene temperature records?, Earth System Dynamics 7, 419 (2016).
  * Fredriksen and Rypdal [2017] H. B. Fredriksen and M. Rypdal, Long-range persistence in global surface temperatures explained by linear multibox energy balance models, Journal of Climate 10.1175/JCLI-D-16-0877.1 (2017).
  * Fraedrich  _et al._ [2004] K. Fraedrich, U. Luksch, and R. Blender, 1/f-model for long-time memory of the ocean surface temperature, Physical Review E 70, 037301 (2004).
  * Rypdal  _et al._ [2013] K. Rypdal, L. Østvand, and M. Rypdal, Long-range memory in Earth’s surface temperature on time scales from months to centuries, Journal of Geophysical Research: Atmospheres 118, 7046 (2013).
  * Zhu  _et al._ [2019] F. Zhu, J. Emile-Geay, N. P. McKay, G. J. Hakim, D. Khider, T. R. Ault, E. J. Steig, S. Dee, and J. W. Kirchner, Climate models can correctly simulate the continuum of global-average temperature variability, Proceedings of the National Academy of Sciences of the United States of America 116, 8728 (2019).
  * Mandelbrot and Van Ness [1968] B. B. Mandelbrot and J. W. Van Ness, Fractional Brownian Motions, Fractional Noises and Applications, SIAM Review 10, 10.1137/1010093 (1968).
  * Malamud and Turcotte [1999] B. D. Malamud and D. L. Turcotte, Self-affine time series: measures of weak and strong persistence, J. Stat. Plan. Infer. 80, 173 (1999).
  * Peixto and Oort [1984] J. P. Peixto and A. H. Oort, Physics of climate, Reviews of Modern Physics 56, 365 (1984).
  * Rohling et al. [2012] E. J. Rohling et al., Making sense of palaeoclimate sensitivity, Nature 491, 683 (2012).
  * Rohling  _et al._ [2018] E. J. Rohling, G. Marino, G. L. Foster, P. A. Goodwin, A. S. von der Heydt, and P. Köhler, Comparing Climate Sensitivity, Past and Present, Annual Review of Marine Science 10, 261 (2018).
  * Zhang  _et al._ [2019] F. Zhang, Y. Qiang Sun, L. Magnusson, R. Buizza, S. J. Lin, J. H. Chen, and K. Emanuel, What is the predictability limit of midlatitude weather?, Journal of the Atmospheric Sciences 10.1175/JAS-D-18-0269.1 (2019).
  * Pelletier [2002] J. D. Pelletier, Natural variability of atmospheric temperatures and geomagnetic intensity over a wide range of time scales, Proceedings of the National Academy of Sciences of the United States of America 99, 2546 (2002).
  * Lohmann and Ditlevsen [2018] J. Lohmann and P. D. Ditlevsen, Random and externally controlled occurrences of Dansgaard–Oeschger events, Climate of the Past 14, 609 (2018).
  * Lovejoy and Lambert [2019] S. Lovejoy and F. Lambert, Spiky fluctuations and scaling in high-resolution epica ice core dust fluxes, Climate of the Past 15, 1999 (2019).
  * Marvel  _et al._ [2016] K. Marvel, G. A. Schmidt, R. L. Miller, and L. S. Nazarenko, Implications for climate sensitivity from the response to individual forcings, Nature Climate Change 6, 386 (2016).
  * Schurer  _et al._ [2013] A. P. Schurer, G. C. Hegerl, M. E. Mann, S. F. B. Tett, and S. J. Phipps, Separating Forced from Chaotic Climate Variability over the Past Millennium, Journal of Climate 26, 6954 (2013).
  * Lovejoy and Schertzer [2012] S. Lovejoy and D. Schertzer, Stochastic and scaling climate sensitivities: Solar, volcanic and orbital forcings, Geophysical Research Letters 39, https://doi.org/10.1029/2012GL051871 (2012).
  * Lovejoy and Varotsos [2016] S. Lovejoy and C. Varotsos, Scaling regimes and linear/nonlinear responses of last millennium climate to volcanic and solar forcings, Earth System Dynamics 10.5194/esd-7-133-2016 (2016).
  * Rypdal and Rypdal [2016] K. Rypdal and M. Rypdal, Comment on “Scaling regimes and linear/nonlinear responses of last millennium climate to volcanic and solar forcing” by S. Lovejoy and C. Varotsos (2016), Earth System Dynamics 7, 597 (2016).
  * Mitchell [1976] J. Mitchell, An overview of climatic variability and its causal mechanisms, Quaternary Research 6, 481 (1976).
  * Henderson-Sellers and McGuffie [1987] A. Henderson-Sellers and K. McGuffie,  _A Climate Modelling Primer_ (Wiley, Chichester, UK, 1987).
  * PAGES 2k Consortium. [2019] PAGES 2k Consortium., Consistent multidecadal variability in global temperature reconstructions and simulations over the Common Era, Nature Geoscience 12, 643 (2019).
  * [47] See the Supplemental Material for further information on data sets and benchmarks.
  * Otto-Bliesner  _et al._ [2016] B. L. Otto-Bliesner, E. C. Brady, J. Fasullo, A. Jahn, L. Landrum, S. Stevenson, N. Rosenbloom, A. Mai, and G. Strand, Climate variability and change since 850 CE an ensemble approach with the Community Earth System Model, Bulletin of the American Meteorological Society 10.1175/BAMS-D-14-00233.1 (2016).
  * Jungclaus et al. [2010] J. H. Jungclaus et al., Climate and carbon-cycle variability over the last millennium, Climate of the Past 10.5194/cp-6-723-2010 (2010).
  * Zhong  _et al._ [2018] Y. Zhong, A. Jahn, G. H. Miller, and A. Geirsdottir, Asymmetric Cooling of the Atlantic and Pacific Arctic During the Past Two Millennia: A Dual Observation-Modeling Study, Geophysical Research Letters 45, 12,412 (2018).
  * Braconnot  _et al._ [2019] P. Braconnot, D. Zhu, O. Marti, and J. Servonnat, Strengths and challenges for transient Mid- to Late Holocene simulations with dynamical vegetation, Climate of the Past 15, 997 (2019).
  * Fischer and Jungclaus [2011] N. Fischer and J. H. Jungclaus, Evolution of the seasonal temperature cycle in a transient Holocene simulation: Orbital forcing and sea-ice, Climate of the Past 10.5194/cp-7-1139-2011 (2011).
  * Liu [2009] Z. Liu, Transient simulation of last deglaciation with a new mechanism for Bølling–Allerød warming, Science 325, 10.1126/science.1171041 (2009).
  * Morice  _et al._ [2012] C. P. Morice, J. J. Kennedy, N. A. Rayner, and P. D. Jones, Quantifying uncertainties in global and regional temperature change using an ensemble of observational estimates: The HadCRUT4 data set, Journal of Geophysical Research Atmospheres 117, 10.1029/2011JD017187 (2012).
  * Hersbach  _et al._ [2020] H. Hersbach, B. Bell, P. Berrisford, and S. Hirahara, The ERA5 global reanalysis, Quarterly Journal of the Royal Meteorological Society , 1999 (2020).
  * PAGES2k Consortium. [2017] PAGES2k Consortium., A global multiproxy database for temperature reconstructions of the Common Era, Scientific Data 10.1038/sdata.2017.88 (2017).
  * Kantelhardt [2011] J. W. Kantelhardt, Fractal and multifractal time series, in  _Mathematics of Complexity and Dynamical Systems_, edited by R. A. Meyers (Springer, New York, 2011) pp. 463–487.
  * Laepple  _et al._ [2018] T. Laepple, T. Münch, M. Casado, M. Hoerhold, A. Landais, and S. Kipfstuhl, On the similarity and apparent cycles of isotopic variations in East Antarctic snow pits, Cryosphere 10.5194/tc-12-169-2018 (2018).
  * Casado  _et al._ [2020] M. Casado, T. Münch, and T. Laepple, Climatic information archived in ice cores: impact of intermittency and diffusion on the recorded isotopic signal in Antarctica, Climate of the Past 16, 1581 (2020).
  * Crowley [2000] T. J. Crowley, Causes of climate change over the past 1000 years, Science 10.1126/science.289.5477.270 (2000).
  * Hegerl  _et al._ [2007] G. C. Hegerl, T. J. Crowley, M. Allen, W. T. Hyde, H. N. Pollack, J. Smerdon, and E. Zorita, Detection of Human Influence on a New, Validated 1500-Year Temperature Reconstruction, Journal of Climate 20, 650 (2007).
  * Braconnot  _et al._ [2012] P. Braconnot, S. P. Harrison, M. Kageyama, P. J. Bartlein, V. Masson-Delmotte, A. Abe-Ouchi, B. Otto-Bliesner, and Y. Zhao, Evaluation of climate models using palaeoclimatic data, Nature Climate Change 2, 417 (2012).
  * Schmidt et al. [2012] G. A. Schmidt et al., Climate forcing reconstructions for use in PMIP simulations of the Last Millennium (v1.1), Geoscientific Model Development 5, 185 (2012).
  * Delaygue and Bard [2011] G. Delaygue and E. Bard, An Antarctic view of Beryllium-10 and solar activity for the past millennium, Climate Dynamics 36, 2201 (2011).
  * Steinhilber  _et al._ [2009] F. Steinhilber, J. Beer, and C. Fröhlich, Total solar irradiance during the Holocene, Geophysical Research Letters 36, 10.1029/2009GL040142 (2009).
  * Wang  _et al._ [2005] Y. Wang, J. L. Lean, and N. R. Sheeley, Jr., Modeling the Sun’s Magnetic Field and Irradiance since 1713, The Astrophysical Journal 10.1086/429689 (2005).
  * Muscheler  _et al._ [2007] R. Muscheler, F. Joos, J. Beer, S. A. Müller, M. Vonmoos, and I. Snowball, Solar activity during the last 1000yr inferred from radionuclide records, Quaternary Science Reviews 26, 82 (2007).
  * Vieira and Solanki [2010] L. E. Vieira and S. K. Solanki, Evolution of the solar magnetic flux on time scales of years to millenia, Astronomy and Astrophysics 10.1051/0004-6361/200913276 (2010).
  * Gao  _et al._ [2008] C. Gao, A. Robock, and C. Ammann, Volcanic forcing of climate over the past 1500 years: An improved ice core-based index for climate models, Journal of Geophysical Research: Atmospheres 113, 10.1029/2008JD010239 (2008).
  * Berger [1978] A. L. Berger, Long-term variations of daily insolation and Quaternary climatic changes., Journal of Atmospheric Sciences 10.1175/1520-0469(1978)035¡2362:ltvodi¿2.0.co;2 (1978).
  * Crucifix [2016] M. Crucifix, Palinsol - Package (R) (2016).
  * Toohey and Sigl [2017] M. Toohey and M. Sigl, Reconstructed volcanic stratospheric sulfur injections and aerosol optical depth, 500 BCE to 1900 CE, version 2, World Data Center for Climate (WDCC) at DKRZ 10.1594/WDCC/eVolv2k_v2 (2017).
  * Fröhlich [2006] C. Fröhlich, Solar irradiance variability since 1978: Revision of the PMOD composite during solar cycle 21, Space Science Reviews 125, 53 (2006).
  * Keeling  _et al._ [1976] C. D. Keeling, R. B. Bacastow, A. E. Bainbridge, C. A. Ekdahl Jr., P. R. Guenther, L. S. Waterman, and J. F. S. Chin, Atmospheric carbon dioxide variations at mauna loa observatory, hawaii, Tellus 28, 538 (1976).
  * Pongratz  _et al._ [2008] J. Pongratz, C. Reick, T. Raddatz, and M. Claussen, A reconstruction of global agricultural areas and land cover for the last millennium, Global Biogeochemical Cycles 22, 10.1029/2007GB003153 (2008).
  * Myhre  _et al._ [1998] G. Myhre, E. J. Highwood, K. P. Shine, and F. Stordal, New estimates of radiative forcing due to well mixed greenhouse gases, Geophysical Research Letters 25, 2715 (1998).
  * Schmidt et al. [2011] G. A. Schmidt et al., Climate forcing reconstructions for use in PMIP simulations of the last millennium (v1.0), Geoscientific Model Development 4, 33 (2011).
  * Myhre et al. [2014] G. Myhre et al., Anthropogenic and natural radiative forcing supplementary material., in  _Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change_, edited by Intergovernmental Panel on Climate Change (Cambridge University Press, Cambridge, 2014) pp. 1–30.
  * Laepple and Huybers [2013] T. Laepple and P. Huybers, Reconciling discrepancies between Uk37 and Mg/Ca reconstructions of Holocene marine temperature variability, Earth and Planetary Science Letters 375, 418 (2013).
  * Chatfield and Xing [2019] C. Chatfield and H. Xing,  _The Analysis of Time Series: An Introduction with R_ , 7th ed., Chapman & Hall / CRC Texts in Statistical Science (CRC Press, 2019).
  * Fredriksen and Rypdal [2016] H. B. Fredriksen and K. Rypdal, Spectral characteristics of instrumental and climate model surface temperatures, Journal of Climate 29, 1253 (2016).
  * Rehfeld  _et al._ [2018] K. Rehfeld, T. Münch, S. L. Ho, and T. Laepple, Global patterns of declining temperature variability from the Last Glacial Maximum to the Holocene, Nature 554, 356 (2018).
  * Percival and Walden [1993] D. B. Percival and A. T. Walden,  _Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques_ (Cambridge University Press, Cambridge, UK, 1993).
  * Yiou  _et al._ [1996] P. Yiou, E. Baert, and M. F. Loutre, Spectral analysis of climate data, Surveys in Geophysics 17, 619 (1996).
  * Kirchner [2005] J. W. Kirchner, Aliasing in $1/f^{a}$ noise spectra: Origins, consequences, and remedies, Phys. Rev. E 71, 66110 (2005).
  * Østvand  _et al._ [2014] L. Østvand, T. Nilsen, K. Rypdal, D. Divine, and M. Rypdal, Long-range memory in internal and forced dynamics of millennium-long climate model simulations, Earth System Dynamics 5, 295 (2014).
  * Fleiss and Cohen [1973] J. L. Fleiss and J. Cohen, The Equivalence of Weighted Kappa and the Intraclass Correlation Coefficient as Measures of Reliability, Educational and Psychological Measurement 33, 613 (1973).
  * Geoffroy  _et al._ [2013] O. Geoffroy, D. Saint-Martin, D. J. Olivié, A. Voldoire, G. Bellon, and S. Tytéca, Transient climate response in a two-layer energy-balance model. Part I: Analytical solution and parameter calibration using CMIP5 AOGCM experiments, Journal of Climate 26, 1841 (2013).
  * MacMynowski  _et al._ [2011] D. G. MacMynowski, H.-J. Shin, and K. Caldeira, The frequency response of temperature and precipitation in a climate model, Geophysical Research Letters 38, https://doi.org/10.1029/2011GL048623 (2011).
  * Kirkevåg  _et al._ [2008] A. Kirkevåg, T. Iversen, J. E. Kristjánsson, Ø. Seland, and J. B. Debernard, On the additivity of climate response to anthropogenic aerosols and CO2, and the enhancement of future global warming by carbonaceous aerosols, Tellus, Series A: Dynamic Meteorology and Oceanography 60 A, 10.1111/j.1600-0870.2008.00308.x (2008).
  * Shiogama  _et al._ [2013] H. Shiogama, D. A. Stone, T. Nagashima, T. Nozawa, and S. Emori, On the linear additivity of climate forcing-response relationships at global and continental scales, International Journal of Climatology 33, 10.1002/joc.3607 (2013).
  * Meehl  _et al._ [2004] G. A. Meehl, W. M. Washington, C. M. Ammann, J. M. Arblaster, T. M. Wigley, and C. Tebaldi, Combinations of natural and anthropogenic forcings in twentieth-century climate, Journal of Climate 17, 10.1175/1520-0442(2004)017¡3721:CONAAF¿2.0.CO;2 (2004).
  * Ramaswamy and Chen [1997] V. Ramaswamy and C. T. Chen, Linear additivity of climate response for combined albedo and greenhouse perturbations, Geophysical Research Letters 24, 10.1029/97GL00248 (1997).
  * Jungclaus [2020] J. Jungclaus, private communication (2020).
  * Danabasoglu [2008] G. Danabasoglu, On Multidecadal Variability of the Atlantic Meridional Overturning Circulation in the Community Climate System Model Version 3, Journal of Climate 21, 5524 (2008).
  * Kunz and Laepple [2021] T. Kunz and T. Laepple, Frequency-dependent estimation of effective spatial degrees of freedom, Journal of Climate 34, 10.1175/JCLI-D-20-0228.1 (2021).
  * Kunz and Laepple [2018] T. Kunz and T. Laepple, Time-scale dependent estimation of spatial degrees of freedom, in  _EGU General Assembly Conference Abstracts_ , EGU General Assembly Conference Abstracts (2018) p. 17995.
  * Yeh  _et al._ [2009] S.-W. Yeh, J.-S. Kug, B. Dewitte, M.-H. Kwon, B. P. Kirtman, and F.-F. Jin, El Niño in a changing climate, Nature 461, 511 (2009).
  * Reschke  _et al._ [2019] M. Reschke, K. Rehfeld, and T. Laepple, Empirical estimate of the signal content of Holocene temperature proxy records, Climate of the Past 15, 521 (2019).
  * Lücke  _et al._ [2019] L. J. Lücke, G. C. Hegerl, A. P. Schurer, and R. Wilson, Effects of memory biases on variability of temperature reconstructions, Journal of Climate 32, 10.1175/JCLI-D-19-0184.1 (2019).
  * Fleiss  _et al._ [2013] J. L. Fleiss, B. Levin, and M. C. Paik,  _Statistical Methods for Rates and Proportions_, 3rd ed., Wiley Series in Probability and Statistics (Wiley, 2013).
  * Clauset  _et al._ [2009] A. Clauset, C. R. Shalizi, and M. E. J. Newman, Power-law distributions in empirical data, SIAM Review 51, 661–703 (2009).
  * Kirtman et al. [2012] B. P. Kirtman et al., Impact of ocean model resolution on CCSM climate simulations, Climate Dynamics 39, 1303 (2012).
  * Klavans  _et al._ [2017] J. M. Klavans, A. Poppick, S. Sun, and E. J. Moyer, The influence of model resolution on temperature variability, Climate Dynamics 48, 3035 (2017).
  * Hodson and Sutton [2012] D. L. Hodson and R. T. Sutton, The impact of resolution on the adjustment and decadal variability of the Atlantic meridional overturning circulation in a coupled climate model, Climate Dynamics 39, 10.1007/s00382-012-1309-0 (2012).
  * Gray et al. [2010] L. J. Gray et al., Solar influences on climate, Reviews of Geophysics 48, 10.1029/2009RG000282 (2010).
  * Hertwig  _et al._ [2015] E. Hertwig, F. Lunkeit, and K. Fraedrich, Low-frequency climate variability of an aquaplanet, Theoretical and Applied Climatology 121, 459 (2015).
  * Peltier and Vettoretti [2014] W. R. Peltier and G. Vettoretti, Dansgaard-Oeschger oscillations predicted in a comprehensive model of glacial climate: A ”kicked” salt oscillator in the Atlantic, Geophysical Research Letters 41, 7306 (2014).
  * Klockmann  _et al._ [2020] M. Klockmann, U. Mikolajewicz, H. Kleppin, and J. Marotzke, Coupling of the Subpolar Gyre and the Overturning Circulation During Abrupt Glacial Climate Transitions, Geophysical Research Letters 47, 10.1029/2020GL090361 (2020).
  * Fernández-Donado  _et al._ [2013] L. Fernández-Donado, J. González-Rouco, C. Raible, C. Ammann, D. Barriopedro, E. García-Bustamante, J. H. Jungclaus, S. Lorenz, J. Luterbacher, S. J. Phipps, _et al._ , Large-scale temperature response to external forcing in simulations and reconstructions of the last millennium, Climate of the Past 9, 393 (2013).
  * Kunz  _et al._ [2020] T. Kunz, A. M. Dolman, and T. Laepple, A spectral approach to estimating the timescale-dependent uncertainty of paleoclimate records – Part 1: Theoretical concept, Climate of the Past 16, 1469 (2020).
  * Dolman  _et al._ [2021] A. M. Dolman, T. Kunz, J. Groeneveld, and T. Laepple, A spectral approach to estimating the timescale-dependent uncertainty of paleoclimate records – part 2: Application and interpretation, Climate of the Past 17, 825 (2021).
  * [113] https://github.com/paleovar/timescaledependency.
  * Fournier [2021] A. Fournier, private communication (2021).
  * Mallat [2009] S. Mallat, Chapter 6 - wavelet zoom, in  _A Wavelet Tour of Signal Processing_ (Academic Press, Boston, 2009) 3rd ed., pp. 205–261.
  * Stratonovich [1967] R. L. Stratonovich,  _Topics in the Theory of Random Noise_ (Gordon and Breach, New York, 1967).
  * Lindner [2006] B. Lindner, Superposition of many independent spike trains is generally not a poisson process, Phys. Rev. E 73, 022901 (2006).
  * Papale [2018] P. Papale, Global time-size distribution of volcanic eruptions on Earth, Scientific Reports 8, 6838 (2018).