The Structure of Climate Variability Across Scales

1.1 Scales in the Climate System

One of the fascinating aspects of the climate system is the close relationship between the spatial and temporal scales of the relevant physical processes. This accounts for the success of scaling analyses of the equations of motion and the systematic derivation of simplified versions of the primitive equations, such as the quasi-geostrophic or the shallow-water equations (e.g., Franzke et al., 2019; Klein, 2010; Majda & Wang, 2006; Vallis, 2017). For instance, the quasi-geostrophic equations are valid in the limit of a small Rossby number (Vallis, 2017) and describe Rossby and synoptic-scale waves and, thus, provide an excellent conceptual model to understand many important aspects of the atmosphere and ocean.

The many physical processes in the Earth's climate system span a vast dynamic range, both in space (from 10⁻³ to 10⁷ m) and time (from seconds to millions of years) (Figure 2). Williams et al. (2017) provide a census of atmospheric processes, the variability of which range from seconds to decades. In the climate system, we typically deal with the following physical processes and associated scales: turbulent eddies on time scales of a few seconds and length scales of millimeters to centimeters, convective activity on temporal scales of hours and spatial scales of hundreds of meters to a few kilometers, synoptic weather systems varying diurnally on spatial scales of hundreds to thousands of kilometers, large-scale teleconnection patterns with an intraseasonal to interannual temporal variability and spatial scales that can span an entire hemisphere, the coupled atmosphere-ocean system which varies from decadal to centennial time scales and a global spatial scale, and the ice ages that represent global variations on millennial time scales (Figure 2). The main four components of the climate system (atmosphere, ocean, land, and cryosphere) tend to operate on different time scales that interact nonlinearly with each other creating a plethora of interesting effects and feedbacks (Peters et al., 2004; Rial et al., 2004; Williams et al., 2017).

An intriguing property of the climate system is that despite the fact that we have to deal with many different physical processes, the variability constitutes a continuum of fluctuations, that is, while the variability spectrum may be interspersed by spikes belonging to some particular and well-defined forcing process (e.g., daily, annual, or Milankovich cycles), the vast part of the spectrum is continuous and scales over large ranges.

1.2 Power Law Scaling

By scaling, we mean the power law relationship between the amplitude of fluctuations and their probability of occurrence on a given temporal or spatial scale:

(1)

where f is an arbitrary function which can either be deterministic or stochastic, y is a climate variable or time, and γ denotes the scaling exponent, a factor which allows us to zoom in and out. In case of f being a stochastic function, the equality has to be interpreted as equality in distribution. When considering a time series, f is a stochastic process and equation 1 implies that the variability of short time scales is statistically similar to the variability on longer time scales. This also implies that no preferred time scale exists. Furthermore, this equation describes a self-similar process (Lamperti, 1962); if y would denote time, then equation 1 would imply that the variance would go to infinity for increasing time scales. Furthermore, the fact that climate data exhibit scaling indicates that the statistical properties remain independent of the scale (Hurst, 1951; Feder, 1988; Franzke et al., 2012; Kolmogorov, 1940; Lamperti, 1962; Mandelbrot & Van Ness, 1968; Mandelbrot, 1982; Taqqu, 2013) as is the case for fractals (Feder, 1988). The scaling property might already be a familiar concept from power spectrum analyses where, in addition to pronounced peaks, one also examines for the existence of linear slopes in a double logarithmic scale representation (e.g., Huybers & Curry, 2006; Wunsch, 2003).

In Figure 3 we display time series sample paths in order to illustrate the scaling property; these were generated from an Autoregressive Fractional Integrated Moving Average (ARFIMA) scaling model (see section 12 and Appendix A). The displayed long-range dependence process has a slope of 0.75 in a log-log plot of fluctuation function versus scale, while a short-range dependence (SRD) process has a slope of 0.5 at long time scales. A slope of 0.5 corresponds to white noise which means that the process is uncorrelated (Figure 3c). The power spectrum (Figures 3d) exhibits the corresponding behavior of increasing power for lower frequencies (with a singularity at zero) of a long-range dependence process exhibiting while the SRD spectrum is flat at low frequencies. Scaling in intensities is displayed in Figure 4 for the α-stable distribution.

1.3 Climate Variability Across Scales

The first attempt to conceptualize atmospheric variability over a wide range of scales has been made by Mitchell (1976). Mitchell's ambitious composite spectrum (Figure 5) ranged from hours to the age of the Earth and focused on the peaks in the power spectrum, thus emphasizing the quasiperiodic phenomena in the climate system and its forcings. Although Mitchell (1976) made a candid admission that his spectrum was mostly an “educated guess,” and despite subsequent improvements in climate and paleoclimate data, the original work has achieved almost iconic status.

Mitchell's scale-bound view led to a climate dynamics framework that emphasizes the importance of numerous processes occurring at well-defined time scales and the separation into quasiperiodic “foreground” processes (illustrated as sharp peaks in Figure 5) and the “unimportant background noise.” We argue that while this division is not wrong per se, it can only explain a small fraction of the overall variability and the underlying climate system dynamics. Wunsch (2003) showed that the quasiperiodic signals represent only a small fraction of the total variability which is more akin to a Lorentzian spectrum of an autoregressive process, while Pelletier (1997) and Huybers and Curry (2006) put an emphasis on the power law behavior of the background spectrum.

See Figure 2a of Lovejoy (2015b) for an illustration of the scaling regimes.

Some recent studies focused more on the continuum aspects of the spectra (Huybers & Curry, 2006; Pelletier, 1998; Paillard, 2001). For instance, Huybers and Curry (2006) reported qualitatively similar results for the macroweather and climate regimes, while Nilsen et al. (2016) provided quantitative evidence that supports the hypothesis of just one scaling regime at least for the Holocene. Nilsen et al. (2016) also question whether it is meaningful to classify climate variability into universal regimes on time scales where we observe forced global climate changes and in particular geological time scales. The reason is that the variability on the long time scales is fundamentally forced by time-dependent external processes, for example, the Milankovich cycle; hence, its statistics are time varying (Nilsen et al., 2016). On shorter temporal scales, on the other hand, scaling is better established in many climatic data sets for a wide range of spatial and intensity ranges. Furthermore, it has been recognized that quasiperiodic signals represent only a small fraction of the total climate variability, and while many studies have focused on understanding these quasiperiodic signals, we argue that the continuous variance spectrum is of equal significance and deserving of future research efforts.

1.4 Scope of the Review

Because scales and scaling properties in the climate system are hard to adequately cover in a single paper, we will restrict this review to topics relevant to the interpretation and reconstruction of time series and to the impacts of scaling on climate variability, trends, prediction, and climate sensitivity. While we cover the potential physical mechanisms behind scaling, we can only provide a broad and nonrigorous introduction to the mathematical framework of scaling processes. More rigorous treatments can be found elsewhere (Beran, 1994; Beran et al., 2013; Baillie, 1996; Doukhan et al., 2002; Embrechts & Maejima, 2007; Guegan, 2005; Lovejoy & Schertzer, 2013; Palma, 2007; Samorodnitsky, 2007, 2016).

Our review is structured as follows: Section 6 covers the basic ideas of scaling and estimation methods; section 26 provides empirical evidence of scaling in climatic time series; section 33 discusses applications of scaling-like trend detection, climate prediction, and climate sensitivity. We end with an outlook and open research questions in section 37.

2.1 Scaling and Power Laws

2.1.1 Scaling From Dimensional Analysis

In the physical sciences, scaling is a well-known and long established concept (Bolster et al., 2011; Longair, 2003; Watkins et al., 2016). For instance, scaling can be used to explain: (i) how a pendulum's angular frequency depends on its length or (ii) how the gravitational force between two bodies depends on their distance from one another.

In the first example, the angular frequency ω depends on the length l as

(2)

where g is the gravitational acceleration. In the second example, Newton's law of universal gravitation states that the gravitational force, F, between two bodies with masses m₁ and m₂, is inversely proportional to the square of the distance between their centers, r, as

(3)

where G is the universal gravitational constant. The scaling property in both examples is so well established that it can be used to extrapolate and to test the behavior of systems outside their initial observable range. It can easily be seen that equations 2 and 3 are different forms of the power law from equation 1 with γ equal to −(1/2) and −2, respectively.

While originally a result of empirical observation, the above equations can also be derived from dimensional analyses. This embodies the physical principle of similarity, which requires that (natural) physical laws should be independent of (human) physical units used to describe a system. According to Buckingham's Π theorem (Buckingham, 1914; Meinsma, 2019), dimensional analysis can be used to show that any physical equation involving n variables can be rewritten using n-m dimensionless parameters, where m ≥ 0, thus revealing possible scaling relations which can then be empirically tested (Bolster et al., 2011).

Dimensional analysis remains a very powerful technique for systems which resist analytic or numerical treatment. The prime example is geophysical fluid turbulence. In 1941 Kolmogorov (Kolmogorov, 1991b, 1991a) derived a scaling relationship between turbulent kinetic energy E and the horizontal scale as measured by wavenumber k for isotropic turbulence. Thereby, he derived the Kolmogorov -5/3 spectrum(for details and underlying assumptions see Vallis, e.g., 2017):

(4)

While a power law distribution of the energy spectrum has been confirmed by observational evidence in the atmosphere (Nastrom & Gage, 1985; Straus & Ditlevsen, 1999), the exact exponent is still a matter of debate (Lovejoy et al., 2007; Lovejoy & Schertzer, 2013). For instance, Lovejoy et al. (2007) have shown that the atmosphere is anisotropic with different scaling exponents in the horizontal and vertical directions, which violates Kolmogorov's assumption of isotropy. Also, the theoretical −(5/3) scaling for large horizontal scales is −2.4 according to aircraft measurements (Lovejoy et al., 2009). This does not invalidate dimensional analysis but only shows that some of the underlying assumptions made by Kolmogorov in his first model (homogeneous and isotropic three-dimensional turbulence) describe an idealized system but are typically not valid in the real atmosphere or ocean, where vertical stratification, jet streams, and the presence of boundaries prevents full isotropy and homogeneity.

Another example of scaling is the addition of N random numbers, where the standard error scales as σ_N∝N^1/2, a result familiar to all scientists from the undergraduate laboratory and the treatment of experimental errors (e.g., Wilks, 2011). Interestingly, this result can be connected to a physical situation, by considering the root-mean-square of the displacement y_N from the origin of the first N steps of a random walk, which is one of the most basic stochastic models for a time series. In a typical one-dimensional discrete random walk, a particle may start at a location and each step moves it either to the left or to the right with equal probability. The resulting root-mean-square of the total displacement, y, after N steps scales with N^1/2 which can also be expressed in terms of time, t, as

(5)

This describes the growth of the diffusing edge of a particle cloud executing Brownian motion (See Appendix B) (Bouchaud & Potters, 2003). The random walk model is statistically self-similar; that is, the time series generated by a random walk looks approximately the same as parts of it. In other words, the shapes and behaviors of the time series are independent of the time scale under consideration. Mathematically, statistical self-similarity can be written as

(6)

and is equivalent to the scaling relationship in equation 1 where refers to that both sides are equally distributed. Here, γ_SS is the self-similarity parameter. In some processes, such as fractional Brownian motion, this is identical to the Hurst exponent H. The Hurst exponent H is named after Harold Edwin Hurst who first identified a scaling relationship investigating the flow levels of the Nile river and other reservoirs (Doukhan et al., 2002; Hurst, 1951, 1957). He developed the R/S method (see details below in Appendices 41 and D) to estimate the scaling exponent. A list of used exponents is given in Table 1.

Table 1. Table of Scaling Exponents

Exponent	Name	Relationship to other exponents
γ	general power law exponent
γSS	self-similarity exponent
H	Hurst exponent	where H measures long-
		range dependence
α	stability exponent
β	power spectrum exponent from a station-	β:=2H−1 where H measures
	ary process	long-range dependence
d	long-range dependence parameter	for Gaussian processes
τ(q)	multifractal exponent/Renyi scaling
	exponent

• Note. d is used in the statistics community in autoregressive fractional integrated moving average models. These models are asymptotically self-similar. H is used in the physical and climatological communities and can be a measure of long-range dependence or self-similarity in systems with Gaussian fluctuations. Here, we use H only as a measure of long-range dependence.

The range of problems we can handle with scaling analysis can be greatly broadened if we introduce the concept of fractals by considering scaling exponents γ which are nonrational. Just as in the integer or rational cases, there is physically instructive information in fractal exponents that can go beyond that from dimensional analysis. These nonrational exponents will play an important role from now on since they are necessary to describe the observed scaling in climate time series due to long-range dependence and heavy-tailed probability density functions (PDFs). They will be discussed in the following subsections.

2.2 Scaling in PDFs and Non-Gaussianity

2.3 Long-Range Dependence

Long-range dependence is characterized by a slow, power law decay of the autocorrelation function. This implies that even long ago states still affect the current state, thus, even far apart in time states, show dependence on each other.

The most basic long-range dependence model is the fractional Brownian motion (See Appendix B). The main difference between fractional Brownian motion and regular Brownian motion is that in Brownian motion the increments are independent of each other while in fractional Brownian motion such increments are dependent in time (Figure 6). This dependence actually covers the whole past; that is the reason that this model is sometimes also called long-term persistence or long memory (for H>0.5). There are different definitions of fractional Brownian motion, and we refer to the specialist literature for more details (e.g., Beran et al., 2013; Beran, 1994; Embrechts & Maejima, 2007; Lévy, 1953; Mandelbrot & Van Ness, 1968).

While fractional Brownian motion is a continuous-time process, the statistics literature prefers a more flexible model, the discrete time ARFIMA (e.g., Beran, 1994; Hosking, 1981; Granger, 1978; Granger & Joyeux, 1980):

(7)

where B denotes the back shift operator BX_t=X_t−1,B²X_t=X_t−2,…. The polynomials Φ and Ψ are defined as and , where p and q are integers and denote the order of the autoregressive Φ and moving average Ψ parts, respectively. The noise variables Z_t are assumed to be independent Gaussian distributed with zero mean and constant variance . See Appendix A for more details.

However, the ARFIMA model can also be generalized to use α-stable distributed increments (Franzke et al., 2012; Graves et al., 2017; Kokoszka & Taqqu, 1994; Stoev & Taqqu, 2005). For these infinite variance models no agreed upon definition of long-range dependence exists (Samorodnitsky, 2016). Note that for d=0 the ARFIMA model reduces to the Autoregressive Moving Average model which is a SRD process. In general, ARFIMA models can also be driven by non-Gaussian (e.g., t-distributed) noise (Graves et al., 2017). ARFIMA models are more flexible than fractional Brownian motion since they combine a long-range dependence component with SRD behavior (Beran, 1994; Beran et al., 2013; Franzke et al., 2012; Graves et al., 2015). The R package ARFIMA can be used to estimate ARFIMA models (Veenstra, 2012).

These are the two most important and widely used paradigmatic models of long-range dependence, but since they were not derived from basic physical laws their use in climate research was originally, and continous to be, met with criticism (e.g., Klemes, 1974; Maraun et al., 2004; Mann, 2011). Long-range dependence also implies that even the most distant past still influences the current and future climate, which appears at odds with common intuition. Many geophysical equations of motion such as the Navier-Stokes or the primitive equations are usually Markovian, that is, their current state only depends on the immediately preceding state and not on states in the more distant past. Furthermore, they do not have memory terms (Chorin & Hald, 2013; Gottwald et al., 2017; Mori, 1965; Zwanzig, 1973, 2001). This fact appears to be at odds with the observed (non-Markovian) long-range dependence behavior of many climate time series and has led to much debate (Bunde et al., 2014; Cohn & Lins, 2005; Franzke, 2012; Maraun et al., 2004; Mann, 2011; Percival et al., 2001; Vyushin & Kushner, 2009). The debate stems from the fact that the underlying equations of motion are Markovian. However, long-range dependence is frequently seen in time series from an aggregated system rather than data from a less ambiguous physical variable, and so the apparent paradox may be illusory since even Markovian systems can appear non-Markovian when not observing the full system. We will discuss possible physical mechanisms to explain this behavior in section 13.

2.4 Multifractals

In section 1, we discussed scaling in precipitation intensities and in temperature time series. For intensity fields as well as time series, there are notions of multifractality that generalize self-similar scaling.

Intensity fields in geophysics can have spatial characteristics that are consistent with random cascades (Kahane & Peyriere, 1976; Kantelhardt, 2009; Sornette, 2006). In such cascades, the intensity in a spatial region distributes nonuniformly between its smaller-scale subregions according to multiplicative processes. The simplest example is the binomial cascade introduced by Kahane (1985). This model originates in turbulence theory, as a rigorous analysis of the Kolmogorov-Obukhov model for spatial variability of the energy dissipation rate (Kolmogorov, 1962; Obukhov, 1962). The multiplicative chaos model (Riedi et al., 1999) is a modern version of the same idea.

The binomial cascade and the multiplicative chaos models define singular (nonsmooth) measures. By construction, the qth moments of the region-averaged intensities are power laws in spatial scale, with exponents that depend concavely on q. Consequently, the distributions of intensities between different spatial regions become increasingly leptokurtic with decreasing scale.

A multifractal time series X(t) is one where the qth moment of an increment |X(t+Δt)−X(t)| scales with the time lag Δt, with an exponent ζ(q) that depends concavely on q. The scaling function ζ(q) is linear for self-similar processes.

There are several ways to construct multifractal stochastic processes from multifractal measures. In most constructions, a multifractal intensity field on the time axis determines the amplitudes in the time series, analogous to how the energy dissipation rate determines the amplitude of velocity field fluctuations in turbulence theory.

For strictly concave scaling functions the distributions of increments are more leptokurtic on short time scales than on longer time scales. Consequently, all multifractal time series are non-Gaussian. The reverse implication does not hold. It is well known that unless one carefully verifies scaling of higher-order moments, standard techniques for estimation of multifractality can lead to spurious results for time series with non-Gaussian marginal distributions.

While multifractals are an abstract concept, they are useful for modeling time series with volatility clustering in time series, where the serial correlations between large and small amplitude events are different.

Applications of multifractal models in climate science have been shown by Schmitt et al. (1995). More recently, Ashkenazy et al. (2003) analyzed climate data from the past 100 kyr and found evidence for nonlinearity and clustering of the magnitude of climatic changes, consistent with multifractality. Similar results have been found by Maslov (2014). Evidence of multifractal scaling in temperature, wind, and precipitation has been found by Baranowski et al. (2015), Gan et al. (2007), and Royer et al. (2008). See Appendix E for multifractal estimators.

2.5 Physical Scaling Mechanisms

Scaling, and particularly long-range dependence, is an actively discussed topic in climate research. There is no obvious physical mechanism in the climate system that would allow the distant past to directly affect the current state of the system. Since the equations of motion used in climate models are all usually Markovian and do not contain memory terms, how can we explain the presence of long-range dependence, and scaling, in the climate system?

2.5.1 Model Reduction

Long-range dependence can be explained using the Mori-Zwanzig formalism from statistical physics (Gottwald et al., 2017; Mori, 1965; Zwanzig, 1973, 2001) which rigorously demonstrates how model reduction leads to the emergence of memory terms in the reduced equations of motion. Let us consider the following example (Gottwald et al., 2017; Zwanzig, 2001):

(8)

(9)

where L_ij are constant parameters. If we are now only interested in the dynamics of x, we can formally solve for y

(10)

which we can now insert into equation 8

(11)

The first term is a Markovian term from the original equations, the second term is a memory term since it integrates over the past, and the last term is the initial condition which can be considered to be random. This example explicitly shows how one gets memory terms when looking only at parts of the full state vector. Equation 11 is still exactly equivalent to the original system.

Most of our measurements are point measurements or just measurements of a subset of the continuous fields. In either case, their dynamics stem from a low-dimensional system embedded in a climate system of infinite dimensions. The Mori-Zwanzig formalism shows that memory effects arise if only a small part of the full system is observed. Thus, long-range dependence could be a direct result of this observation. While the memory term in equation 11 is fairly general—which makes it impossible to know how exactly memory decays—a power law decay is a possibility, especially when making additional assumptions about the memory kernel. Kupferman (2004) approximated the memory kernel with a power law.

2.5.2 Nonlinearity and Regimes

Lorenz put forward the idea that deterministic systems can be almost intransitive; that is, they can exhibit long-lasting climate changes and hence no unique climate state exists (Lorenz, 1968, 1976). Such long-term anomalies can be a form of scaling in that the variance increases with increasing time scale. Several studies have shown that nonlinearity can lead to scaling (Franzke et al., 2015; Lorenz, 1976; Mesa et al., 2012). Atmospheric circulation regime behavior, a main component of the climate system (Feldstein & Franzke, 2017; Franzke, 2013; Franzke et al., 2011; Hannachi et al., 2017; Ghil & Robertson, 2002; Nicolis, 1990), has been suggested as a prime candidate for scaling (Franzke et al., 2015). An example of atmospheric circulation regimes is given by the quasi-stationary circulation systems like blocking events, which are quasi-stationary high-pressure systems that can last for weeks and cause heat waves and cold spells (Feldstein & Franzke, 2017; Hannachi et al., 2017). It has been shown for very long but finite time series that regime behavior is a plausible mechanism for scaling because the residence times of the regimes are power law distributed (Diebold & Inoue, 2001; Franzke et al., 2015). The residence time is the time the system stays in one regime state. If these time intervals are power law distributed, then the system can exhibit long-range dependence. This implies that memory effects in the climate system may not be needed to explain the apparent scaling of variance with time scale. The origin of this scaling has been found to be associated with the coarse graining of the dynamics into a finite number of specific regimes, leading to non-Markovian dynamics (Nicolis, 1990; Nicolis & Nicolis, 1988, 1995; Nicolis et al., 1997; Vannitsem, 2001).

Recent model experiments suggest also another possible nonlinear mechanism that could explain long-range dependence: the coupling of the atmosphere with other components of the climate system that have very different characteristic time scales. A case in point is ocean-atmosphere coupling, for which a reduced order nonlinear coupled model has been developed recently (De Cruz et al., 2016; Vannitsem et al., 2015). This model employs the quasi-geostrophic equations to describe the large-scale dynamics of the atmosphere and oceans in extratropical regions. The coupling is achieved via an energy balance scheme and momentum transfer through wind stress.

Multiple scaling regimes were found (Figure 7) using a Haar wavelet analysis (see Appendix 44). Remarkably, no low-frequency variability was found in the coupled model for small friction coefficients and the moments peak at a scale of about 10 days and decrease for larger periods. By low-frequency variability, we mean a set of long-periodic, attracting orbits that couple the dynamical modes of the ocean and the atmosphere in this model. If low-frequency variability develops in the system, then additional peaks emerge at 10,000 and 40,000 days. Similar to Lovejoy (2015b), this allows us to define different regimes based on the respective scaling exponents. The structure of the low-frequency variability and long-range dependence critically depends on the water depth (Figures 7b and 7c). This suggests that one plausible explanation of observed scaling regimes lies in the coupling of climate subcomponents. We will further discuss this coupling mechanism in a linear framework next.

2.5.4 Non-Gaussianity and Multiplicative Noise

As discussed above, scaling can also arise from the distribution of the increments or the driving noise in a stochastic process. So far, we only discussed scaling in additive noise processes which in addition may have heavy tails. Also, Gaussian noise can produce power law PDFs when it occurs in a multiplicative or state-dependent process (Bódai & Franzke, 2017; Franzke, 2017; Majda et al., 2009; Penland & Sardeshmukh, 2012; Sornette, 2006; Sardeshmukh & Sura, 2009; Sura & Hannachi, 2015). The simplest multiplicative noise process is the Kesten process (Sornette, 2006), a first-order autoregressive process model with random coefficients:

(12)

where a_n and b_n are independent random variables. Under certain conditions, the Kesten process has a process cumulative probability density function with a power law decay of its tails, that is,

(13)

where γ is the power law exponent.

Stochastic climate theory predicts the presence of multiplicative noise in nonlinear systems (Franzke et al., 2015, 2019, 2005; Franzke & Majda, 2006; Franzke, 2017; Gottwald et al., 2017; Majda et al., 1999, 2001, 2008, 2009; Penland & Sardeshmukh, 2012; Sardeshmukh & Sura, 2009; Sura & Hannachi, 2015). It can also be shown that multiplicative noise leads to power laws over some ranges in stochastic climate models (Majda et al., 2009; Sardeshmukh & Sura, 2009; Sura & Hannachi, 2015). Unlike power law processes, stochastic climate theory also provides mechanisms to limit extremes. This power law roll-off is due to the same nonlinear interaction that causes the multiplicative noise in the first place: the nonlinear interaction between slow and fast components (Franzke et al., 2005; Franzke & Majda, 2006; Majda et al., 1999, 2001, 2008, 2009; Sardeshmukh & Sura, 2009; Sura & Sardeshmukh, 2008). This can be understood as follows: The fast components of the flow, for example, convection, synoptic-scale weather systems, are effectively serially uncorrelated on the time scale of the slow components, for example, Rossby waves or the ocean. This time scale separation allows us to treat the fast components effectively as a noise variable. While there are nonlinear interactions between the slow and the fast components in the climate system, this can be written now as a product of the slow flow variable and a noise variable, that is, multiplicative noise, also called state-dependent noise since the impact of the noise can be modulated by the state of the slow variable. This is consistent with the findings of Sardeshmukh and Sura (2009) where they found evidence in global circulation model simulations that multiplicative noise is due to turbulent adiabatic fluxes and not rapid diabatic forcing fluctuations. An example is wind gusts: If the large-scale wind speed is low, then there are only weak wind gusts; on the other hand, if the large-scale wind speed is high, also, the wind gusts are strong. This behavior can be easily represented by a multiplicative noise where the wind gusts are computed by the product of the large-scale wind speed and a noise. The relevance of multiplicative noise has been shown for sea surface temperature variability (Sura & Sardeshmukh, 2008), atmospheric vorticity variability (Sardeshmukh & Sura, 2009), teleconnection patterns such as the North Atlantic Oscillation (Majda et al., 2009; Önskog et al., 2019), and extreme events (Franzke, 2017; Penland & Sardeshmukh, 2012; Sura, 2013).

While theoretical considerations predict a power law, for example, the generalized central limit theorem (Sornette, 2006), our climate system is of finite size and thus infinitely large events cannot occur which means that the power laws need to cut or roll-off at some intensity or spatial size. This is also consistent with the dynamical systems theory of extremes (Lucarini et al., 2016) which shows that pure power law dynamics cannot occur at arbitrarily large intensities or sizes.

2.6 Estimation Methods for Scaling Exponents

A multitude of estimators have been developed over the years to provide accurate estimates of the scaling exponents, and different estimators infer different aspects of the scaling properties. For instance, most estimators infer the long-range dependence parameter d or the Hurst exponent H of a time series and are insensitive to non-Gaussianity of its amplitudes, which can cause them to differ from the self-similarity parameter γ_SS. When deciding whether or not to use a particular estimator, one should always be aware of the underlying assumptions that went into its construction.

2.7 Spectral Analysis and Its Limitations

Analyzing data by Fourier decomposition is attractive because each sinusoid has an exact and unambiguous time scale: its period. Furthermore, one has to assume that the signal is statistically stationary in time or homogeneous in space for Fourier analysis to be valid. This implies that the probability laws that govern the behavior of the time series do not change over time. As a result, the phases of the Fourier modes are randomly distributed and can be neglected. However, this only applies if the time series is Gaussian because the phases are correlated otherwise; a fact that is often overlooked. If observationally data really had to meet these criteria fully, then Fourier analysis would not be as useful as it has been, so we can relax these conditions in practice to some extend.

A potential problem with Fourier analysis is that the interpretation of the resulting spectra is nontrivial. To illustrate this problem, we consider surrogate time series which have the same multifractal characteristics as the EPICA dust series (Lambert et al., 2008; Lovejoy, 2018). Since the time series are generated by a multifractal stochastic process, even the strong peaks are randomly excited (Figure 8). However, the presence of such strong peaks can be misinterpreted as signatures of important climatic processes (Lovejoy, 2018). This example should also be taken as a cautionary tale against the simple interpretation of peaks in Fourier analysis as quasi-oscillations resulting from physical mechanisms. The important point we want to emphasize here is not that these methods do not work but rather that their resulting spectra need to be carefully interpreted and the model representing the underlying process carefully chosen.

3.1 Station Temperatures

Significant evidence for long-range dependence in station temperature have been reported in many studies (Bunde et al., 2014; Capparelli et al., 2013; Fraedrich & Blender, 2003; Franzke, 2010, 2012; Gao & Franzke, 2017; Graves et al., 2015; Koscielny-Bunde et al., 1998; Ludescher et al., 2015; Løvsletten & Rypdal, 2016, 2015). The exponent H=0.65 of the fluctuation function F(τ)∼τ^H which was determined by DFA was suggested to be universal for local temperature variations. This exponent is related to the exponent β=0.3 in a power spectrum, by β:=2H−1. Figure 9 shows the results for the fluctuation function obtained by DFA (see Appendix 43) for 12 meteorological stations distributed worldwide. In the interannual range (up to the total duration) the stations reveal linear slopes given by an exponent H≈0.65. It is tempting to conclude that atmospheric long-range dependence could be characterized by a single universal value. Later on, it was found that the value H=0.65 might be a transition phenomenon between memory-less inner continents with H=0.5 and some oceanic areas with H≈1 (Fraedrich & Blender, 2003; Ludescher et al., 2015; Løvsletten & Rypdal, 2016).

3.2 Near-Surface Temperatures on Land and Ocean

Observed near-surface temperatures are available since 1900 with a sufficient density to estimate a global pattern of long-range dependence (Jones, 1994; Jones et al., 2001; Parker et al., 1995). Figure 10 shows the power law exponents H for the monthly data estimated by DFA with quadratic trend (DFA2). The results are concentrated in coastal regions, Europe, North America, and along-ship routes. In inner continental locations, long-range dependence is negligible with H≈0.5. Along the coasts the memory increases to H≈[0.6,0.8] and in the central North Atlantic and in the equatorial Indian ocean, the highest values of H≈0.9 are found. This pattern shows that long-range dependence in the considered time range is a marine phenomenon. On land no memory is evident on these time scales far from the coasts. The finding of H≈0.65 in the stations can be explained by the locations of the stations considered in Figure 9. Inhabited areas with observational stations are traditionally along coasts. Note that there can be a huge gradient in long-range dependence as seen along the western North Pacific coast. Later on, Fraedrich and Blender (2003) found that this surface temperature long-range dependence can be found in coupled atmosphere-ocean simulations with a dynamical ocean model but not with a slab ocean model.

A spectrum with β=1 is denoted as 1/f noise (see Appendix F). If this extends to vanishing frequency (or infinite time), stationarity is violated. Clearly, in observational data this limit cannot be attained, and whether the sea surface temperature data are stationary remains open. However, this type of variability indicates that short-term averages, as are typically used in climate science, are not defined and should be interpreted with caution.

The 1/f spectrum of the sea surface temperature in some oceanic regions can be obtained by a diffusion model in a vertical column with two compartments (Fraedrich et al., 2004). A decisive parameter is the diffusivity in the abyssal ocean. This parameter determines turbulent mixing and is caused by tides and the orography. Although the atmospheric forcing is white, the sea surface temperature shows long-range dependence close to nonstationarity. Furthermore, wind also provides a significant part of the mechanical energy for diapycnal mixing in the ocean (Munk & Wunsch, 1998; Wunsch & Ferrari, 2004) which could also lead to this behavior as discussed in section 14.

In summary, the analysis of observational data reveals that long-range dependence in the lower atmosphere is predominantly found in oceanic regions where the variability is close to nonstationarity (1/f spectrum). Far from coasts long-range dependence on decadal time scales is not observed. In transition zones along coasts a spectrum is found which is approximately given by .

3.3 Scaling in Regional and Global Mean Temperatures

Local and regional surface temperature variations have a much greater magnitude and show a weaker scaling than global average temperature variability (e.g., Laepple & Huybers, 2014a). This difference in variability is strongest for interannual and shorter time scales and decreases on longer time scales. The reduced variability of the global mean temperatures reflects cancelation of variability in the global mean, and the weaker cancelation toward lower frequencies is consistent with findings that temperature anomalies have greater spatial autocorrelation toward longer time scales (Jones et al., 1997). This behavior is also reproduced in diffusive energy balance models (Rypdal et al., 2015) where the predicted slope of the spectrum of global mean temperatures is around double the slope of regional temperatures over a large frequency range following from the horizontal diffusive coupling. Thus, it is important to consider the spatial scale analyzed when interpreting the variability and long-term memory of temperature time-series.

3.4 Paleo Data and Simulations

The isotope fraction δ¹⁸O in ice cores can be used as a proxy for the surface temperature due to the different weights of the molecules (Barlow et al., 1997; Dansgaard, 1964). As the annual snowfall is preserved on the Greenland and Antarctic ice sheets, this allows in principle to analyze climate scaling from interannual to multimillennial time scales. However, especially on the shorter time scales, nontemperature effects such as snow redistribution (Fisher et al., 1985), diffusion (Johnsen et al., 2000) and aliasing effects from intermittent snowfall (Laepple et al., 2017) considerably affect the recorded variability and have to be corrected for, to infer about climate scaling (Münch & Laepple, 2018). The oxygen fractions measured in the Greenland ice cores GRIP and GISP2 during the last 10,000 years reveal estimates for H (Blender et al., 2006) which are clearly above 0.5 in the millennial time range (Figure 11). The corresponding exponents in the power spectrum correspond to H≈0.5. Hence, scaling can be assumed, at least approximately. It is remarkable that this result can be obtained in an extremely long coupled atmosphere-ocean simulation (Blender et al., 2006) which reveals intense long-range dependence south of Greenland (see Figure 11) with exponents of similar magnitude but much less long range dependence in other oceanic regions; the Pacific ocean reveals no long-range dependence of a comparable intensity. The simulation was performed under present-day conditions; hence, no external variability is necessary to explain this result. The long-range dependence is related to the variability of the zonally averaged stream function in the North Atlantic Ocean. Evidence for long-range dependence has also been found in other coupled climate models (Fredriksen & Rypdal, 2016; Østvand et al., 2014; Vyushin et al., 2004; Zhu et al., 2010). Vyushin et al. (2004) showed that the inclusion of volcanic eruptions improves the simulation of long-range dependence in climate models.

An important question is how well climate models reproduce the true internal variability on time scales of centuries and longer. The local model surface temperature spectra seem to indicate a lower (or even zero) spectral exponent β for frequencies below 1×10⁻² year⁻¹, but on long time scales the finite sample size errors are so large that this cannot be concluded with high statistical confidence (Fredriksen & Rypdal, 2016). This flattening of the spectra on time scales longer than a century cannot be detected in the instrumental temperatures, since the time series are too short, but they are also not detected in temperature reconstructions of Holocene climate (Laepple & Huybers, 2014b) even after correcting for nonclimate effects on the spectra. On the contrary, some authors claim higher exponents for local temperatures (β≈2.2) for these longer time scales based on composite spectra established from different proxy records (Huybers & Curry, 2006; Lovejoy & Schertzer, 2012b). This conclusion, however, can only be drawn with confidence from proxy records spanning the last glacial period and including the last deglacation and hence may not be valid for the Holocene climate.

Analyzing global mean temperatures that are dominated by external forcing on the other hand suggests that the scaling and variability are comparable between reconstructions and model simulations (Crowley, 2000; Zhu et al., 2019) or climate model simulations may even overestimate global mean temperature variability.

The issue is illustrated in Figure 12a where we have plotted time series of reconstructed annual temperatures for the Northern Hemisphere and the corresponding series derived from the NorESM model with historical forcing. In Figure 12b we have estimated their Haar fluctuation functions (see Appendix 44). The fluctuation level is more than twice as high for the model temperatures on time scales less than 100 years. It turns out that this is due to the higher short-time responses to large volcanic eruptions in the models. This can be seen by elimination of these spikes from the model signal, as shown in the zoomed in signals in Figure 12. The fluctuation function of this “chopped” signal is very close to that of the reconstructed temperature (Figure 12). The cloud of thin curves in Figure 12 is fluctuation functions estimated for 20 realizations of fractional Gaussian noises with Hurst exponent H=0.9 of the same length as the NorESM model run. The width of the cloud suggests that the reconstruction and the chopped model signal are consistent with such a fractional Gaussian noise, although the power at time scales longer than a few centuries is somewhat high. It is easy to verify that this increased power is due to the temperature difference between the Medieval Warming Anomaly and the Little Ice Age. Some authors interpret the power in this oscillation as a signature of a transition to scaling with an exponent β>1 on time scales longer than a few centuries (Huybers & Curry, 2006; Lovejoy & Schertzer, 2012b). On the other hand, Nilsen et al. (2016) argue that existing temperature reconstructions for the Holocene are generally consistent with a single scaling regime with H≈0.9 on all time scales shorter than the duration of the interglacial period.

Similar ideas are advocated by Rypdal and Rypdal (2016b), who demonstrate that temperatures derived from ice cores over the last 100 kyr can be described as sudden transitions between stadials and interstadials superposed on a 1/f noise (H≈1) background. According to this view, monofractal, near Gaussian, scaling with H≈1 is a useful description of the climate noise background in Quaternary climate. Whether a scaling description is appropriate for the succession of transitions between stadials and interstadials is another story and needs more research.

3.5 Scaling Behaviors in Other Meteorological and Climatological Variables

Besides the evidence of long-range dependence in temperatures, there are also scaling behaviors detected in many other variables, including precipitation, river runoff, total ozone, relative humidity, and sea level change. For in situ precipitation records, small long-range dependence parameters have been found by several studies (Jiang et al., 2017; Kantelhardt et al., 2006; Yang, Fu 2019). On average, the DFA exponent H mainly ranges between 0.5 and 0.55, indicating weak long-range dependence. For river runoff much stronger long-range dependence has been detected. Kantelhardt et al. (2006) found that the mean long-range dependence parameter for river runoff is 0.72 based on 42 river runoff records observed from Europe, North and South America, Africa, Australia, and Asia. Wang et al. (2008) detected long-range dependence close to 1/f (H≈1) for the intra-annual Yangtze discharge. As for relative humidity, Chen et al. (2007) reported that the mean DFA exponent H for in situ relative humidity records over China is around 0.75. Recently, there were also results reported indicating that sea level changes are characterized by long-range dependence. The DFA exponent H has a large variation range from 0.60 to 0.95, depending on different regions (Dangendorf et al., 2014). Other variables such as wind speed, atmospheric circulation indices, and ozone anomalies have also been shown to have the scaling behavior (Feng et al., 2009; Franzke et al., 2015; Vyushin et al., 2007; Vyushin & Kushner, 2009; Varotsos & Kirk-Davidoff, 2006).

3.6 Evidence of Multifractal Behavior

Besides long-range dependence that only needs one exponent H to describe monofractal behavior, there is also empirical evidence of multifractal behavior. For instance, in precipitation records, although the measured long-range dependence is weak, pronounced multifractality has been found (Kantelhardt et al., 2006), indicating that precipitation records of different amplitudes have different scaling behaviors. Similar properties also exist in river runoff data (Koscielny-Bunde et al., 2006), where the multifractality was found to be even stronger than that in precipitation records (Kantelhardt et al., 2006). For temperature related records such as the surface mean air temperature and diurnal temperature range, different multifractal behaviors were found over different regions (Lin & Fu, 2008; Yuan, Fu, & Mao, 2013). Such as in the south of the Yangtze River, pronounced multifractality was found in diurnal temperature range records, while in the north, the multifractal behavior is very weak or even nonexistent (Yuan, Fu, & Mao, 2013). Other variables such as wind speed and relative humidity have also been shown to have the multifractal behavior (Baranowski et al., 2015; Kavasseri & Nagarajan, 2005).

4.1 Scaling for Trend Detection

The identification of trends is one of the most frequent and prominent goals in the analysis of geophysical time series (Chandler & Scott, 2011; Wu et al., 2007). Although apparently an easily understandable objective, trend assessment is very challenging, starting with the lack of a precise definition of trend itself. Implicit in the intuitive notion of trend are concepts such as long term, smoothness, or monotonicity, but it is not unambiguously defined how long is “long term” or how smooth needs a pattern to be in order to be a trend. Furthermore, time series characterized by scaling behavior often exhibit features that can be classified broadly as a trend, even in the absence of any genuine trend.

Unlike the notion of trend, stationarity is a well-defined statistical property. A time series (X_t) is weakly stationary if its first and second moments are time invariant (i.e., the mean and variance are constant and the covariance depends only on the time lag between the observations). The trend in a time series can be ascribed to a nonstationary generating process (at least the mean is not constant in time) and described by a trend model. Trend models can be broadly classified as either being (i) deterministic or (ii) stochastic. Deterministic trend models represent deterministic (nonrandom) nonstationary processes which are described by a function evolving in time; one example is trends which are forced by external factors such as anthropogenic greenhouse gas emissions. Stochastic trend models represent stochastic nonstationary processes described by models such as a random walk or an autoregressive integrated moving average model. Such models exhibit apparent trends without any external forcing; instead, these trends are caused by the internal dynamics of the process. Unfortunately, it is extremely difficult, for example, by visual judgment alone (Percival & Rothrock, 2005), to select a trend model or even to distinguish between deterministic and stochastic trend models, particularly in the case of short time series (Franzke, 2012).

Furthermore, weakly stationary processes can generate a time series with an apparent trend, particularly when a short segment of the process is observed. Such “spurious” trends can be misleadingly taken as evidence of nonstationary behavior when in fact the process is stationary. A classical example since long recognized as a potential culprit in the interpretation of climate variability (Wunsch, 1999) is the typical “red noise” (Red noise sometimes means that the spectral power increases with period scale. However, in climate science red noise typically denotes the power spectrum of a first-order autoregressive process which has first increasing power for increasing period but then becomes white noise; also called Lorentzian spectrum. See Figure 3 for an example.) structure of climate records. A red noise (described by a simple first-order autoregressive process) can produce visually appealing trends despite being a stationary process, particularly in the case of short time series.

Long-range dependence processes, for example, described by ARFIMA models, are another type of stationary processes that can produce apparent nonstationary behavior and local “spurious” trends. Since long-range dependence is a feature common in many geophysical time series, a crucial challenge in the identification and estimation of trends is to discriminate between nonstationary processes and stationary long-range dependent processes. The problem is, however, quite difficult since genuine trends generated by a nonstationary process and spurious trends from a long-range dependence process can coexist in the same time series. Disentangling the different contributions to the observed temporal structure is not possible by visual inspection, and even specific methodologies addressing the issue have to rely on substantial assumptions and simplifications, for example, on the type of nonstationary behavior or the dominance of one specific type of process. For instance, the approach proposed by Beran and Feng (2002) of semiparametric fractional autoregressive (SEMIFAR) models considers a trend function modeled nonparametrically, with the remaining components of the model estimated by maximum likelihood. Despite the flexibility of SEMIFAR models, the performance is poor in the case of short time series, and the trend is estimated based on a subjective concept of smoothness. More importantly, discrimination between stochastic and deterministic trends remains difficult to achieve, given that a significant amount of spurious trends associated with long-range dependence behavior alone can be easily included in the nonparametric trend estimation. The statistical test of Berkes et al. (2006) aims to discriminate between stationary long-range dependence time series and nonstationary time series with change points in the mean but requires previous identification of a small number of change points in the time series.

Shifting from a general trend model to a linear trend significantly constrains the problem of trend identification in the presence of long-range dependence. Although such an assumption is hardly realistic and is theoretically limiting, it is nevertheless of practical relevance since the overwhelming majority of trends in geophysical records are reported as the slope from a linear regression model. Most studies are based on the up-front assumption that a time series can be described by a nonstationary linear trend with a stochastic long-range dependence component (e.g., Bunde et al., 2014; Capparelli et al., 2013; Franzke, 2010, 2012; Lennartz & Bunde, 2009; Ludescher et al., 2015; Myrvoll-Nilsen et al., 2019; Rybski & Bunde, 2009) and focus on the assessment of the corresponding uncertainty (e.g., Cohn & Lins, 2005; Koutsoyiannis, 2006; and Koutsoyiannis & Montanari, 2007).

A complementary approach is to test the assumption of a linear deterministic trend itself. The test devised by Phillips and Perron (1988), the PP test, is a classical unit root test for testing nonstationarity in the form of a random walk, which is a scaling process. The test devised by Kwiatkowski et al. (1992), the KPSS test, is a parametric statistical test which assumes as a null hypothesis a deterministic linear trend plus a stationary stochastic noise. Although the two tests can be applied independently, their joint use is recommended for trend assessment purposes (Carrion-i-Silvestre et al., 2001; Kwiatkowski et al., 1992). For a time series with a random walk stochastic trend the PP test should not reject the null hypothesis and the KPSS test should reject the linear trend null hypothesis. Conversely, for a time series with a linear deterministic trend, the PP test should reject the null hypothesis but not the KPSS test. In the case of a times series for which both tests fail to reject the null hypothesis, then the time series or the tests are not sufficiently informative to distinguish between a stochastic (random walk) trend and a deterministic trend. However, if both tests reject their respective null hypothesis, this is an indication that alternative parametrizations for long-term behavior need to be considered, such as long-range dependence. Since long-range dependence is a common feature of geophysical time series, this outcome of rejection of both PP and KPSS test is quite common, for example, in the case of air temperature (Fatichi et al., 2009) or global sea surface temperature (Barbosa, 2011). Unfortunately, these tests require long time series (in order to meet asymptotic assumptions) and are known to have low explanatory power particularly against long-range dependence alternatives (Lee & Schmidt, 1996; Leybourne & Newbold, 1999).

A widely used trend test is the Mann-Kendall test (Kendall, 1948; Mann, 1945), which in its original form is only valid for independent data. The Mann-Kendall test is a nonparametric test which tests for the presence of a monotonic trend without making any assumptions about the form of the trend. This is in contrast to most other trend tests which have to assume some parametric trend form, for example, a linear trend. The Mann-Kendall test has been extended to also account for serial correlation in time series (Hamed & Rao, 1998) and also for the presence of long-range dependence (Hamed, 2008).

A further trend significance test has been developed by Lennartz and Bunde (2009, 2011) and Tamazian et al. (2015). This method has been developed for the DFA method in the presence of long-range dependence in the time series (Bunde et al., 2014; Ludescher et al., 2015). Based on Monte Carlo simulations, they studied how the trend uncertainties vary with the strength of long-range dependence, as well as the data length. This method has been applied to evaluate trend significances of the surface air temperature and the sea ice extent in Antarctica (Bunde et al., 2014; Ludescher et al., 2015; Yuan et al., 2017).

Recent developments in temperature trend significance testing with long-range dependent noise show how we can also incorporate information about forced global temperature changes in the trend estimate (Myrvoll-Nilsen et al., 2019). In that way, one avoids attributing forced changes deviating from, for example, a linear trend as part of the stochastic variability (Gil-Alana, 2005; Fatichi et al., 2009; Franzke, 2012, 2014). Results show that the observed trends since 1900 are significant relative to the noise for most locations and to a larger degree than when assuming a linear trend (Løvsletten & Rypdal, 2016).

4.2 Scaling for Climate Response and Sensitivity

Linear response models, which predict how the climate system will react to a change in forcing, for example, anthropogenic greenhouse gas emissions, have shown considerable success in describing the global temperature response in climate model data, instrumental data, and in multiproxy reconstructions (Caldeira & Myhrvold, 2013; Fredriksen & Rypdal, 2016; Held et al., 2010; Geoffroy et al., 2013; Lovejoy et al., 2015; Østvand et al., 2014; Rypdal & Rypdal, 2014; Rypdal et al., 2015). In particular, Rypdal and Rypdal (2014) demonstrated that a scaling linear response function provides a good description of the global temperature response to radiative forcing over both the historical period and to a multiproxy reconstruction of the temperature over the last millennium.

It has been known for several decades that Coupled Atmosphere-Ocean General Circulation Models exhibit climate responses on multiple time scales (Caldeira & Myhrvold, 2013; Fredriksen & Rypdal, 2017; Held et al., 2010; Geoffroy et al., 2013); that is, there is more than one time constant involved in the response. The scaling response studied in Rypdal and Rypdal (2014) could be considered an approximation to the multiple time scale response, bridging the responses on time scales from years to centuries, and Fredriksen and Rypdal (2017) demonstrates how an energy balance box model can provide such an approximation.

In addition to describing the response to historical radiative forcing, the same scaling response to a white noise stochastic forcing is also consistent with the observed internal variability. One way of extracting the internal variability from observed global temperature is to compute the deterministic, historically forced variability from the response model and subtract this from the observed record. The power spectrum of this estimated internal variability compares well with a power law (Rypdal & Rypdal, 2014). This tendency for a multibox model to form a power law spectrum is studied systematically in Fredriksen and Rypdal (2017) and reflects a well-known result which states that a scaling spectrum can be obtained from the aggregation over an ensemble of first-order autoregressive processes (Granger, 1980); see also section 15. Thus, long-range dependence can be caused by the constructive superposition of SRD processes.

The emergent scale invariance makes it possible to infer equilibrium climate sensitivity from a scaling frequency-dependent climate sensitivity R(f)∼f^β/2. This scaling response implies infinite magnitude response as f→0, and therefore, there must exist a lower-frequency limit of where the scaling response is valid and where the response stabilizes as we go to even lower frequencies. R(f) can be estimated for a given climate model by exploiting the relation between the historic radiative forcing applied to a model and the observed instrumental global temperature. Rypdal et al. (2018) applied this to an ensemble of Earth system models, where the inferred values of R(f) evaluated at f=1/1,000 year⁻¹ correlate strongly to estimates of equilibrium climate sensitivity from idealized model runs. They could use the distribution of estimated R(f) over the model ensemble to constrain the distribution of equilibrium climate sensitivity obtained from the ensemble of idealized runs. Thus, scale-invariant linear response models are useful tools for the estimation of equilibrium climate sensitivity from observation data. The advantage over multibox energy balance models is that the scale-invariant models contain fewer free parameters and are less prone to statistical overfitting.

4.3 Scaling for Climate Prediction

The property of long-range dependence in the climate system raises the question whether models which explicitly include long-range dependence can be used for skillful predictions. The first attempt for climate predictions was tried by Baillie and Chung (2002). Recently, a model was developed for seasonal to decadal predictions, the Stochastic Seasonal to Inter-Annual Prediction System (StocSIPS) (Lovejoy, 2015a; Lovejoy et al., 2015, 2018). StocSIPS is based on the low-frequency limit of a fractional differential equation, the fractional energy balance model. This model is valid for periods between 20 days and 50 years, where intermittency is relatively weak so that a quasi-Gaussian approximation can be used. StocSIPS forecast skill compares favorable with operational long-range forecasting models based on traditional climate models. One advantage of StocSIPS is that data assimilation of observations is not necessary, since it can directly be fitted to observed data. This also implies that downscaling of forecasts is not needed.

Yuan, Fu, and Liu (2013) and Yuan et al. (2014) developed a method for the extraction of the long-range dependence using a fractional integrated statistical model. They proposed a new variable memory kernel which clearly shows how the states from the distant past maintain their impacts over time till the current time. Accordingly, climate variables with long-range dependence can be decomposed into two parts: (i) the memory part, which represents the influences accumulated from the past, and (ii) the residual part, which is related to the current dynamical forcing conditions. With the memory part extracted, one can at least determine on what basis the considered time series will continue to change. By combining this with the estimated residual part, it is possible to make predictions. Therefore, they proposed a new perspective for climate prediction for climate variables with long-range dependence. Because the influence from the past can be extracted quantitatively, one only need to focus on the prediction of the residual part.

Also, statistical models with non-Gaussian features have recently been developed. For instance, Önskog et al. (2018) show that the forecast skill of the North Atlantic Oscillation (Feldstein & Franzke, 2017) increases in a SRD statistical model when non-Gaussian noise is used. Graves et al. (2017) developed a Bayesian framework for ARFIMA models with various non-Gaussian noises and demonstrated its usefulness using the t-stable and the α-stable distributions.