This paper studies the fluctuations of an equilibrium solution reached by a dynamical system under constant external forcing conditions. Applying this definition to climate, the fluctuations of our concern correspond to internal climate variability (ICV) resulting from internal processes without involving any time-varying external forcing. The existence of ICV has been known for a long time (see e.g. a review by Mitchell (1976)). A vast literature has been devoted to the description and the qualification of ICV, and to the understanding of mechanisms responsible for ICV. ICV has different appearances on different temporal and spatial scales. With respect to temporal scales, there was a great excitement generated by the discovery of James and James (1989) that internal dynamics of the atmosphere produce not only variations on the synoptic time scales of days but also variations on much longer time scales of about a few decades. Since then decadal and multi-decadal climate variations have drawn attention of many researchers. With respect to spatial and temporal scales, a variety of modes, ranging from atmospheric teleconnections (Panagiotopoulos et al., (2002) to the coupled modes, such as the El Niño Southern Oscilation (ENSO) (Rasmusson & Carpenter, 1982) in the tropical Pacific, and to the large-scale modes of extra-tropical SST, such as the Pacific Decadal Oscillation (PDO) (Newman et al., 2016) and the Atlantic Multidecadal Oscillation (AMO) (Kerr, 2000), has been investigated using both observations and model simulations.
Despite all the efforts, our understanding about ICV is still incomplete. An indication of the incompleteness can be found in a study on two atmospheric modes simulated by a coupled atmosphere-ocean general circulation model (GCM) (J.-S. von Storch, 1999a). These modes vary in phase with the two parts of the global atmospheric axial angular momentum M on time scales of months to a century. Because of the relation to M, and because M is completely determined by its torque F via ${\scriptscriptstyle \frac{dM}{dt}}=F$ , it was thought that the cause of the low-frequency variations of the two modes can be understood by examining the torque F. This turns out to be a wishful thinking. Against the expectation, the coherence and the phase spectra between M and F break down on time scales longer than a few months (bottom panel in Figure 1), inhibiting establishing any link between M and F on these scales. The result is irritating, since the angular momentum equation ${\scriptscriptstyle \frac{dM}{dt}}=F$ is one of those first principles in which our deterministic thinking is rooted.
In this paper, we aim to clarify the mystery shown in Figure 1 and identify the mechanism responsible for ICV in a possibly general manner. Broadly speaking, there are, to my knowledge, two avenues that can lead to a more general explanation of climate variability. The first one emphasizes the role of nonlinear dynamics and is built on the dynamical systems approach. As reviewed by Dijkstra and Ghil (2005) for low-frequency variability of the ocean circulation, the observed ocean variability can be interpreted as emerging from the road of successive bifurcations leading through multiple equilibria to oscillatory and eventually chaotic solutions. What is less clear with this dynamical systems approach is the extent to which it can be applied to purely internally generated variability. A bifurcation behavior has to be initiated by varying a parameter of the considered system. For the ocean, the parameters studied in this context are wind stress and buoyancy fluxes at the sea surface. As both parameters are also the external forcing of the ocean, the distinction between internally and externally generated variability becomes blurred.
The second avenue emphasizes the stochastic nature and is culminated in the concept of stochastic climate models proposed by Hasselmann (1976). Following this concept, variations of a slow climate component are the response to continuous random excitation by fast climate components (such as the weather). The evolution equation of a slow climate component includes a stochastic forcing representing the fast components and takes the form of a stochastic differential equation. Such evolution equations suggest that the internal climate variability is stochastic and does not result from deterministic responses to a temporally varying external forcing. Despite the wide acceptance of the concept, it is still not definitely clarified what stochasticity means. The meaning of stochasticity is of fundamental importance: A physical phenomenon, such as ICV, cannot be governed by two different types of evolution laws, one formulated in terms of deterministic differential equations and the other in terms of stochastic differential equations. If we accept that all climate components, no matter fast or slow, result from deterministic laws (such as the conservation of angular momentum) that do not contain any stochastic element, a stochastic model can only be an approximation at most. To “upgrade” a stochastic model to something having the same status as the governing deterministic laws, one needs first to show that there is a physical “thing” that escapes the control of the governing deterministic evolution laws.
This papers shows that for dynamical systems of our concern such a “thing” exists. A dynamical system of our concern has the following three properties. First, the system described by the state vector x is governed by the Newtonian type of dynamics such that the temporal evolution of a component x of x is determined by
${\scriptscriptstyle \frac{d(\cdot )}{dt}}$ is the time derivative operator. The set of Eq.(1), formulated for all components of x, represents the deterministic laws governing the evolution of x. f = f(x) is a deterministic non-linear function of x, which contains a constant external forcing as well as a dissipation. A dissipation is needed to counteract the constant external forcing, so that the system reaches an equilibrium state under the constant external forcing, rather than runs away. Hereafter f is referred to as the differential forcing of x, also later when discrete version of Eq.(1) is considered. Eq.(1) represents the classical Newtonian dynamics in the sense that the time rate of change of a “momentum”, i.e. x in Eq.(1), equals the applied “force”, i.e. f(x) in Eq.(1).
Secondly, the system’s state vector x must be multidimensional so that the differential forcing of x, f(x), depends not only on x but also on other components of x. For such a multidimensional system, a differential forcing f of component x at time t describes the interaction of x with other components of x at time t.
Finally, when forced with a constant external forcing for a while, the system reaches an equilibrium state described by a time-varying equilibrium solution. An equilibrium solution fluctuates stationarily, non-periodically, and for ever into infinite time (when left alone). The fluctuations of an equilibrium solution are referred to as the equilibrium fluctuations. Equilibrium fluctuations have finite variances, referred to as the equilibrium variances. Equilibrium variances are determined by the constant external forcing, independent of the initial state, from which the system starts its “journey” toward the equilibrium state. For any component x(t) of an equilibrium solution x(t), x(t) at two time points become increasingly less well correlated with increasing separation between the two time points. As a result, all covariance functions of variables, that are functions of an equilibrium solution x(t) (e.g. the auto-covariance function of a component x of x(t), the auto-covariance function of the differential forcing f of x, and the cross-covariance function between x and f), decay fast in magnitude with increasing time lag.
An example of systems having the just stated properties is an atmospheric (or an oceanic or a coupled atmosphere-ocean) GCM. In this case, the dimension of x is determined by the grid used to discretize the governing equations of motions, which are of the type of Newtonian dynamics and can be transformed into a set of equations in form of Eq.(1). A component x of x is a prognostic variable at one grid point. The differential forcing f of x is a function not only of x but also of other variables at the same and the surrounding grid points. When forcing the GCM with a constant external forcing (e.g. from the sun) for some time, we obtain an equilibrium climate. When rerunning the same GCM under the same external forcing, we obtain the same climate characterized by the same equilibrium variances.
Another example is a gas as a classical many-particle system. Obviously, the governing dynamics are Newtonian and can be described by a set of equations in form of Eq.(1). A component x of x is the position or velocity of one particle. The differential forcing f of x is a function not only of x but also of positions and velocities of other particles. When putting a gas in contact with a heat reservoir as an external forcing, it reaches an equilibrium state characterized by a constant mean kinetic energy of the gas molecules. Independent of the initial states of the individual gas molecules, a gas brought in contact with the same heat reservoir reaches the same equilibrium state characterized by the same mean kinetic energy.
A final example is the 3-dimensional Lorenz model (Lorenz, 1963). The model with the parameters given in Lorenz (1963) is derived from a description of convective motions within a fluid layer. The external forcing, which represents the difference in temperature at the upper and lower boundaries of the layer, is merged into one of the model parameters. The Lorenz model is governed by three equations in form of Eq.(1). The differential forcing f of one component x is a function not only of x but also of the other components of the state vector x. With the same parameter (corresponding to the same external forcing), two equilibrium solutions starting from two different initial states fluctuate in similar manner (top row of Figure 2) and have the same variances. The auto-covariance functions of x (bottom panel in Figure 5) and f (top panel in Figure 6), as well as the cross-covariance function between x and f (bottom panel in Figure 6), decay fast in magnitude with increasing time lag.
To this end we note that the system of our concern is a simpler and idealized version of the reality. Take the climate as an example. The real climate, which is always subjected to some time-varying external forcing (e.g. to a changing atmospheric greenhouse gas concentration that leads to a changing radiative forcing of the system), is generally not in an equilibrium state. A climate subjected to a constant forcing is an artificially constructed and simpler version of the real climate. We believe that the real climate cannot be fully understood without first fulling understand this idealized climate.
We note further that even though some properties of equilibrium fluctuations resemble those produced by a chaotic dynamical system, we intend not to equal the systems of our concern with chaotic dynamical systems (Eckmann & Ruelle, 1985). The main reason is that the latter has been an acronym for properties associated with an attractor, which are not the concern of the present study.
In this paper we should use the Lorenz model to illustrate and support our ideas. We do so not only because the Lorenz model belongs to the dynamical systems of our concern, but also because of the limited computer resources. Given that the number of simulated years per day of a typical modern climate model is around 50 (Balaji et al., 2017), and given that we do need at least multicentury-long, if not millennia-long, simulations for studying the climate variability at the lowest frequencies, it is not feasible to produce a sufficiently large number of sufficiently long solutions of a climate model with the present computer facility. Such excessive computations, which are needed to avoid sampling errors, such as those found in Figure 1, can be carried out for the 3-dimensional Lorenz model.
The paper is organized as follows. After some preliminary notions introduced in Section 2, Section 3 discusses two requisites needed for the considerations in Section 4. The first one is a clean definition of spectrum (Section 3.1). The second is a specification of the low-frequency shape of spectrum for variables that are functions of an equilibrium solution x(t) (Section 3.2). We show that for these variables, whose covariance functions decay fast in magnitude with increasing time lag, the spectra of x and f and the cross-spectrum between x and f are white at the lowest frequencies. In Section 4, equipped with the spectrum definition identified in Section 3.1, we derive the known spectral relations between x and f arising from the differential operator in Eq.(1). Equipped with the low-frequency spectral shape identified in Section 3.2, we are able to point out that for component x having a white low-frequency spectral plateau, the variance related to the spectral plateau, that makes up a non-zero portion of the total equilibrium variance of x, cannot be determined by f, since f de facto ceases to fluctuate at near-zero frequencies. Section 5 discusses ideas suggesting why x fluctuates at the lowest frequencies, and what, if not the differential forcing f, determines the low-frequency spectral plateau of x. Conclusions and a discussion on the implication of this paper is provided in Section 6.
Since a system described by Eq.(1) can generally only be solved numerically, we consider a discretized version of Eq.(1),
where s ∈ ℤ is an integer and t = sδt indicates the time at the s-th time step. For the sake of simplicity, we set the integration time step δt to one and consider,
f is also known as the first difference of x. The limit δt → 0 turns out to be irrelevant and will be briefly discussed at a later point. Hereafter, for a clear distinction between discrete and continuous variables, we denote a discrete variable by a subscript, e.g. x_{s}, and a continuous variable by an argument, e.g. x(t). For an easy indexing, we write a finite discrete record of a variable x as
which contains total 2N+1 time steps. Here and hereafter, the subscript _{N} denotes a quantity related to a finite record. In theory, N in Eq.(4) can go to infinity, since an equilibrium solution can, when left alone, be extended to infinity. In practice, N is always finite.
Two operators are needed. The first one is
where E_{m} is an average over quantities related to m equilibrium solutions. The m equilibrium solutions correspond to m points in the respective phase space. They are obtained by integrating the dynamical system from m initial states that are different but belong to the same equilibrium state. Thus, E(∙) takes all reachable equilibrium solutions of the considered equilibrium state into account. Note that E(∙) is purely deterministic. The link of E(∙) to the expectation operator ɛ(·) used in probability theory is discussed in Appendix A.
The second operator is the Fourier transform $\mathcal{F}\{\cdot \}$ , or the discrete Fourier transform when dealing with discrete variables. The Fourier transform $\mathcal{F}\{\cdot \}$ of a discrete real-valued function α indexed by integer $\tau \in \mathbb{Z},\mathcal{F}\left\{\alpha \right\}$ , is defined as
which operates on infinite values of α^{τ} and results in a function of ω defined on the real interval [–1/2, 1/2] (H. von Storch & Zwiers, 1999). If α is a solution of a dynamical system, τ indexes the time step. If α is a covariance function, τ indexes the time lag. $\mathcal{F}\left\{\alpha \right\}$ is generally a complex-valued function, but becomes real-valued, if ατ is symmetric about τ = 0. As pointed out by Priestley (1981), $\mathcal{F}\left\{\alpha \right\}$ does not always exist. A necessary condition for the existence of $\mathcal{F}\left\{\alpha \right\}$ is that α is absolutely summable, meaning
A variable that is a function of an equilibrium solution is not absolutely summable. The spectrum that decomposes the variance of this variable in frequency domain must be defined as the Fourier transform of the auto-covariance function of this variable.
This section seeks for equilibrium solutions of a system of our concern, a) a general definition of spectrum, and b) a specification of the general low-frequency shape of spectra of variables that are functions of an equilibrium solution.
Throughout this subsection and if not specified otherwise, x denotes a variable related to an equilibrium solution, i.e. x can be a component of an equilibrium solution or a differential forcing. We first concentrate on the definition of auto-spectrum and extend the definition to cross-spectrum at the end of this subsection.
The spectrum of x decomposes the variance of x in frequency domain. In order to concentrate on variance, we consider throughout this paper only anomaly variables obtained by subtracting the respective mean from the variable. We assume, and as can be expected from many Newtonian type systems, that the variance is ergodic, i.e.
with
x_{s} in Eq.(9a) is taken from a finite record {x_{s,N}} indexed as in Eq.(4), with s running from –N to N. ${x}_{s}^{i}$ in Eq.(9b) indicates the value of x associated with the i-th equilibrium solution at an arbitrary time point s. Since x is a function of an equilibrium solution, the exact value of s in Eq.(9b) is irrelevant. In words, ergodicity in variance means that the variance obtained by averaging ${x}_{s}^{2}$ over a sufficiently large number of time points equals the variance obtained by averaging ${x}_{s}^{2}$ over a sufficiently large number of ensemble points. Figure 3 shows that for the three Lorenz components, the variance is indeed ergodic.
The ergodicity in variance induces one to think that the spectrum is also ergodic, and can be derived by applying the discrete Fourier transform either to an auto-covariance function defined using a time average, or to an auto-covariance function defined using an ensemble average. The former is given by
with
and the latter by
with
The two lines in Eq.(11) differ only with respect to the values taken by the running index s. Appendix B (below Eq.(B14)) explains why 1/(2N + 1 – |τ|) is chosen as the normalization factor in Eq.(11), even though other factors are used in statistics literature (H. von Storch & Zwiers, 1999). As in Eq.(9a), x_{s} in Eq.(11) is taken from a finite record {x_{s,N}}. As in Eq.(9b), the choice of s in Eq.(13) is irrelevant for x that is a function of an equilibrium solution.
Both c_{τ} and γ_{|} result from an average of two-time products x_{s–t}x_{s}. Similar to the variance that can be obtained by averaging the product ${x}_{s}^{2}$ either over time points or over ensemble points, an auto-covariance at lag τ can be obtained by averaging the two-time product x_{s–τ}x_{s} either over time points following Eq.(11) and Eq.(10), or over ensemble points following Eq.(13) and Eq.(12). However, there is a big difference between the two averages. The former is an average over different numbers of two-time products depending on τ, since the number of available two-time products in a single record {x_{s,N}} equals 2N – |τ| + 1 and hence decreases with increasing |τ|. On the contrary, the latter is an average over the same number of two-time products, regardless of the value of τ.
Assuming that auto-covariance function exists, either the limit given in Eq.(10) or the limit given in Eq.(12) must converge. If for any value of τ, the two limits converge to the same value, the auto-covariance function is ergodic. However, c_{τ,N} with large values of τ cannot converge in the limit N → ∞ due to the decrease in the number of available two-time products with increasing time lag. For any value of N in the limit N → ∞, c_{τ,N} at time lags close to N is always averaged over a limited number of two-time products. For τ = 2N, c_{τ,N} reduces to the “raw” two-time product x_{–N}x_{N} without being subjected to any average! Thus, it is the lack of available two-time products at large time lags that makes it impossible for c_{τ,N} to converge to zero with increasing N, even when c_{τ,N} would converge to zero, if there were sufficiently many two-time products available for averaging. To ensure the convergence, we must increase the available two-time products, which is only possible by employing more equilibrium solutions.
For the Lorenz system, Figure 4 shows how c_{τ,N} behaves with increasing N; and Figure 5 shows how γ_{τ,m} behaves with increasing m. Since both c_{τ,N} and γ_{τ,m} are symmetric about τ = 0, only positive lags are shown. At small time lags with τ ≤ 1000 (left columns in Figures 4 and 5), c_{τ,N} converges with increasing N towards γ_{τ,m} obtained with the largest value of m: c_{τ,N} obtained with 2N = 10^{6} (bottom left panel in Figure 4) is visually identical to γ_{τ,m} obtained with m = 10^{6} (bottom left panel in Figure 5). However, we need – for defining a spectrum – the values of c_{τ,N} and γ_{τ,m} at all τ with τ ∈ ℤ, since the Fourier transform is an operation on auto-covariance function at all τ with τ ∈ ℤ (see Eq.(6)). When considering all time lags (right columns in Figure 4 and Figure 5), the behavior of c_{τ,N} with increasing N is clearly different from the behavior of γ_{τ,m} with increasing m. While γ_{τ,m} becomes increasingly smooth and the magnitudes of γ_{τ,m} at large time lags converge to zero with increasing m, c_{τ,N} at large time lags remains erratic without any sign of convergence, reflecting the above described lack of two-time products available for time averaging.
Given that only γ_{τ,m} converges for all τ ∈ ℤ, whereas c_{τ,N} cannot converge for large values of τ ∈ ℤ, we are left with no other choice than defining spectrum as the Fourier transform of γ obtained by averaging over all reachable equilibrium solutions. Another way to see this is to recall that a necessary condition for the Fourier transform of an infinite series to exist is that the series is absolutely summable. c_{τ} is not absolutely summable, since the magnitude of c_{τ} remains large at large time lags due to the lack of two-time products available for averaging. Generally, the spectrum obtained from one equilibrium solution is erratic and differs from the spectrum obtained from another equilibrium solution.
Hereafter, a spectrum, denoted by Γ(ω), will be defined as the Fourier transform of γ_{τ}
Since Γ(ω) is symmetric, Γ(ω) over non-negative frequencies is considered throughout this paper. Appendix B shows that Γ(ω) of a variable x decomposes the variance of x as defined in Eq.(8). Appendix B also shows that Γ(ω) can be equally expressed as
where S_{ω,N} is the sample spectrum (periodogram) defined in terms of Fourier coefficients used to Fourier decompose a finite record {x_{s,N}} (see Eq.(B5)). In this paper, all spectra of the Lorenz system are estimated following Eq.(15), rather than Eq.(14), since Γ(ω) can be better approximated using E_{m}(S_{ω,N}) with a moderate value of m than using auto-covariance function γ_{τ,m} with the same value of m. The situation comes about since the Fourier transform is sensitive to the values of an auto-covariance function at large time lags. Even when an auto-correlation function would decay to zero when sufficiently averaged, an estimate of this auto-correlation function, obtained by averaging over limited number of two-time products, has generally a error structure such that it varies slowly around zero at large time lags (Bartlett, 1946; Bloomfield, 1976). To accurately estimate spectrum as the Fourier transform of an auto-covariance function, this error structure must be sufficiently suppressed, which requires considering a large value of m. Figure 5 shows that γ_{τ,m} is still clearly non-zero at large time lags when employing an ensemble of size m = 10^{4} (second panel in the right column), and only converges visually to zero for m larger than 10^{5} (third and fourth panel in the right column) . In contrast to that, Appendix B (Figure B1) shows for the Lorenz components that E_{m}(S_{ω,N}) converges fairly fast with m and we are able to obtain a stable spectrum estimate already with m = 10^{3}.
To further reinforce our conclusion that spectrum should be defined by taking all equilibrium solutions into account, Appendix C shows for the Lorenz system that the sample spectrum S_{ω,N}, that is defined for a single equilibrium solution, does not converge with increasing N. The convergence of S_{ω,N} in the limit N → ∞ is needed to obtain a spectrum that decomposes the variance as defined in Eq.(9a).
Before leaving this subsection, we extend the definition given in Eq.(14) for an auto-spectrum to the definition of a cross-spectrum between two variables, e.g. the cross-spectrum between a component x and its differential forcing f, Γ^{xf} (ω). Γ^{xf} (ω) is defined in analogy to Eq.(14) as the Fourier transform of the cross-covariance function
Similar to Eq.(15), Γ^{xf} (ω) is identical to the double limit of the sample cross-spectrum ${S}_{\omega ,N}^{xf}$ . The latter is defined in analogy to Eq.(B5) in terms of Fourier coefficient of x multiplied with the complex conjugate of the Fourier coefficient of f.
A spectrum Γ(ω) defined in Eq.(14) has the spectral value at ω = 0
Γ(0) is linked to the low-frequency shape of Γ(ω). In this subsection, Γ(ω) can be the auto-spectrum of a variable (e.g. a component x or its differential forcing f), or the cross-spectrum between two variables (e.g. between x and f). The variables are functions of an equilibrium solution whose properties are specified in Section 1. For such variables, covariance functions γ_{τ} decay fast in magnitude with increasing time lag. This is illustrated for the Lorenz system in the bottom panel of Figure 5 for the auto-covariance functions of x, ${\gamma}_{\tau}^{x}$ , and in Figure 6 for the auto-covariance functions of f, ${\gamma}_{\tau}^{f}$ , and for the cross-covariance functions between x and f, ${\gamma}_{\tau}^{xf}$ .
For a decaying covariance function, there exists a finite time lag τ^{†}, hereafter referred to as the cut-off time lag, and a negligibly small positive ϵ, such that
We consider ϵ being negligibly small when Γ(ω) in Eq.(14) can be well approximated by
Knowing the value of τ^{†} that allows the above approximation, we can define (when considering only posivite frequencies) the frequency interval
with
Ω covers the low-frequency end of the entire non-negative frequency interval [0,1/2]. The exact value of ω_{0}, the right end of Ω, and with that the exact width of Ω, depend on τ^{†}, which in turn depends on the accuracy of Eq.(19). For γ_{τ} whose magnitude converges sufficiently fast to zero, τ^{†} is finite for a negligibly small ϵ. The width of Ω is not zero. The near-zero frequencies ω ϵ Ω are physically meaningful for equilibrium solutions, which can extend to infinite times and fluctuate on infinite timescales.
For ω ϵ Ω and τ <τ^{†}, we have 2πωτ < 2πωτ^{†} ≪ 1, so that e^{–i2πωτ} in Eq.(19) can be well approximated by one, since all higher-order terms in the Taylor expansion of e^{–i2πωτ} are negligibly small relative to one. We obtain
Thus, the low-frequency part of Γ(ω) is essentially white and equal the respective spectral value at frequency zero, Γ(0).
For a spectrum whose low-frequency part is white, Γ(0) not only represents the spectral value at a single frequency, ω = 0, but also determines the white spectral level at the lowest frequencies and the variance related to this white low-frequency spectral level.
With the spectrum definition given in Eq.(14), we are now in the position to translate Eq.(3) formulated in time domain into spectral relations between x and f formulated in frequency domain. Such a relation can take the form of the cross-spectrum between x and f, Γ^{xf}, or a relation between the spectrum of x, Γ^{x}, and the spectrum of f, Γ^{f}. We derive the spectral relations following H. von Storch and Zwiers (1999, chapter 11.4.4). In their case, in which the spectral relation between a stochastic process X_{t} and its first difference X_{t}–X_{t–1} is considered, a covariance function is indisputably defined in terms of an expectation operator ɛ(·) that involves joint probability functions (s. Appendix A). In our case, in which a dynamical system in its equilibrium state is considered, we have the choice of defining a covariance function either as c_{τ} using just one equilibrium solution, or as γ_{τ} using all reachable equilibrium solutions. Once the choice of γ_{τ} is made, for the reasons given in the last section, the derivation can be carried out in the same way as in H. von Storch and Zwiers (1999) by replacing ɛ(·) by E(·).
We obtain, by shifting s in Eq.(3) to s–τ, multiplying the shifted Eq.(3) with Eq.(3), and then applying E(·),
and by multiplying Eq.(3) with x_{s–1–τ}, and applying E(·),
Applying Fourier transform to γ^{x}, γ^{f}, and γ^{fx} in Eq.(23) and Eq.(24) yields
with the proportionality factors
where the following equality is used
Thus, both Γ^{f} (ω) and Γ^{xf} (ω) are proportional to the spectrum of x, Γ^{x} (ω).
Γ^{xf}(ω) is a complex-valued function and can be described by its amplitude, known as the amplitude spectrum A^{xf}(ω), or the coherence and phase spectrum, C^{xf}(ω) and P^{xf}(ω). C^{xf}(ω) is a normalized version of the amplitude of Γ^{xf}(ω), defined as (A^{xf}(ω))^{2}/(Γ^{x}(ω)Γ^{>f}(ω)). P^{xf}(ω) equals the arctangent of the ratio of the imaginary to the real part of Γ^{xf}(ω). Given Γ^{xf}(ω) in Eq.(26), these spectra take the following forms
Since the period of the tangent function is π, P^{xf} (ω) is not unique and given by two expressions in Eq.(30). Since A^{xf}(ω) and C^{xf}(ω) are symmetric and P^{xf}(ω) is anti-symmetric, only positive frequencies need to be considered.
For sufficiently small ω satisfying 2πω ≪ 1, G^{f}(ω) and G^{xf}() can be well approximated by
Eq.(25), Eq.(26), Eq.(29) and Eq.(30) reduce, for 2πω ≪ 1, to
The above spectral relations are known (H. von Storch & Zwiers, 1999; Jenkins & Watts, 1968; J.-S. von Storch, 1999b). They result from the first difference operator (·)_{s} – (·)_{s–1} as a high-pass filter, characterized by the proportionality factor G^{f}(ω) (also known as the gain function). J.-S. von Storch (1999b) discussed the case when Γ^{x}(ω) is white at low frequencies. However, these relations are generally valid for solutions of any dynamical system that is governed by evolution equations in form of Eq.(1), independent whether the solution fluctuates stationary or not, and regardless how covariance functions of the solution behave with increasing time lag. In such a general case, Γ^{x}(ω) can have any spectral shape at the lowest frequencies, including having a power law Γ^{x}(ω) ∝ ω^{–n} with n ≥ 2. Γ^{f}(0) would then be finite and non-zero if n = 2, and infinite if n > 2.
The situation is different when considering equilibrium solutions of a dynamical system of our concern specified in Section 1. In this case, all three covariance functions, ${\gamma}_{\tau}^{x},{\gamma}_{\tau}^{f}$ and ${\gamma}_{\tau}^{xf}$ , have a finite cut-off time lag. Denote the largest of the three by τ^{†}. We can define a low-frequency interval Ω having a non-zero width (as it is done in Eq.(20)), such that for ω ϵ Ω, i.e. for ω ⋍ 0 but not strictly ω = 0, Γ^{x}(ω), Γ^{f}(ω), and Γ^{xf}(ω) are essentially white. Γ^{x}(ω) in Eq.(33)–Eq.(35) can then be replaced by a finite Γ^{x}(0) for ω ϵ Ω, noting that Γ^{x}(0) must be finite since the respective equilibrium variance is finite. Eq.(33)–Eq.(35) can then be replaced by
At frequency ω = 0, we have exactly
Thus, different from Eq.(33) which represents only a relation between Γ^{x}(ω) and Γ^{f}(ω) at the lowest non-zero frequencies, Eq.(37) specifies the spectral values of Γ^{f}(ω) at frequency ω ϵ Ω: f de facto ceases to fluctuate at frequency ω ϵ Ω due to the finite constant spectral level of Γ^{x} (and the smallness of ω), even though for any non-zero ω ϵ Ω, Γ^{f}(ω) is numerically not exactly zero and decreases with decreasing ω at the rate of ω^{2}. Because f de facto ceases to fluctuate at frequency ω ϵ Ω, Γ^{xf}(ω) and A^{xf}(ω) are de facto zero for ω ϵ Ω, even though for any non-zero ω ϵ Ω, Γ^{xf}(ω) is numerically not exactly zero and whose amplitude decreases with decreasing ω at the rate of ω, and the normalized amplitude spectrum, the coherence spectrum C^{xf}(ω), is one at all non-zero frequencies.
That Γ^{f}(ω) is de facto zero has a strong physical implication. As it is shown in Figure 7a for the Lorenz system, the spectrum of component x, Γ^{x}(ω), is non-zero at the lowest frequencies. Figure B2 further shows how Γ^{x}(ω) converges to a non-zero Γ^{x}(0), as ω approaches zero. With Γ^{x}(0) ≠ 0 and Ω having a non-zero width, the spectral plateau of Γ^{x}(ω) over the frequency range Ω makes up a non-zero portion of the total equilibrium variance of x. This portion of the total equilibrium variance of x, hereafter referred to as Γ(0)-variance of x, cannot be generated by fluctuations of the differential forcing f of x, since f de facto ceases to fluctuate at ω ϵ Ω.
The above consideration applies only for solutions with infinite length. In practice where only finite records are available., complications arise. Figure 7 shows for the three Lorenz components, estimates of the spectra of x (top left), the spectra of differential forcing f (top right), and the coherence and phase spectra between x and f (bottom). All three Lorenz components reveal similar spectral behavior. The spectra of f have a maximum at frequencies around several 10^{–2} cycles per time step (cpts), and decrease thereafter with decreasing frequency, with the rate of decrease being close to ω^{2}, consistent with Eq.(37). The decrease of spectrum of f is stopped at frequencies around 10^{–5}–10^{–4} cpts. Thereafter and down to the smallest resolved frequency, the spectra of f are white. The coherence spectrum increases from zero at the highest frequencies, at which the spectral values of f are extremely small, to one at intermediate frequencies, consistent with Eq.(31). The phase at these intermediate frequencies is –90°, consistent with Eq.(30). At frequencies lower than 10^{–5}–10^{–4} cpts, at which the spectra of f stop to decrease with decreasing frequency, the coherence and –90°-phase relation completely break down.
The low-frequency behaviors shown in Figure 7 are closely linked to the frailty of f. Appendix D shows that a spectrum estimated using finite records is subjected to truncation errors. Appendix D shows further that the spectrum of truncation error is white at the lowest frequencies. As the spectrum of f decreases with ω and becomes smaller than the truncation error at sufficiently low frequencies, the truncation error dominates. The spectrum of f approaches the white error spectrum. At these low frequencies, the coherence breaks down and the phase deviates strongly from -90°. For the Lorenz system, the frailty of f at the lowest frequencies can be more clearly seen by plotting the spectra in linear (rather than log) scale. Figure 8 shows Γ^{f}(ω) and A^{xf}(ω) in linear scale. Since Γ^{f}(ω) at the lowest frequencies, at which f is weaker than the white truncation errors, is at least five to six orders of magnitude smaller than Γ^{f}(ω) at frequencies around 10^{–2} cycles per time step (cpts), Γ^{f}(ω) and A^{xf}(ω) at the lowest frequencies appear to be zero relative to Γ^{f}(ω) and A^{xf}(ω) at frequencies around 10^{–2} cpts.
Note that the truncation error is not critical for the spectrum of x at the lowest frequencies. Relative to the high low-frequency spectral plateau of x, the truncation error can always be made negligibly small by considering a not too small N.
Before leaving this section, we point out that the spectral behaviors at the lowest frequencies found for the Lorenz system (Figure 7) resemble those found for the global axial angular momentum M and the global torque F simulated by a coupled GCM (Figure 1). In both cases, albeit that the spectra obtained from the GCM are much noisier, the behaviors of both models are characterized by a) an essentially white low-frequency plateau of spectrum of M or x, b) a tendency of decreasing spectral level of F or f with decreasing frequency before the decrease is stopped at the lowest frequencies by the truncation error, and c) the lack of coherence between M and F or between x and f at the lowest frequencies. The similarity suggests that these spectral features are general properties of a Newtonian type dynamical system subjected to similar truncation errors.
Generally, spectra are derived using limited data in practice. The resulting spectrum estimate is subjected both to sampling error and to truncation error. The former arises when applying the average E_{m} to an ensemble of too small size. When no ensemble is available, one often uses an average over several chunks of a record, or a smoothing in frequency domain. For instance, the spectra in Figure 1 were obtained by cutting the available simulation into 5 pieces. The truncation error on the other hand emerges when estimating near-zero spectral values at the lowest frequencies using records of finite length. The spectra in Figure 1 are contaminated by both types of errors. The spectra shown in Figure 7, which are derived by using ensembles of size 1000, are essentially free of sampling error and mostly only contaminated by the truncation error.
We point out another side of the problem arising from using finite records, often encountered in climate applications. Generally, a prognostic variable describing a climate in an equilibrium state should reveal a low-frequency spectral plateau. Using records with limited length impedes fully resolving such spectral plateaus. For this reason, white low-frequency plateaus of prognostic climate variables (in an equilibrium climate state) are not always perceived.
Finally, we point out that the time derivative operator ${\scriptscriptstyle \frac{d(\cdot )}{dt}}$ filters out low-frequency spectral signals in the same way as the first difference (·)_{s} – (·)_{s}–1. When translating Eq.(1) into a spectral relation in frequency domain, we find that Γ^{f} is proportional to Γ^{x} in a way similar to that given in Eq.(25), but with the proportionality factor G^{f}(ω) = (2sin(πω))^{2} being replaced by (2πω)^{2}, which approaches G^{f}(ω) = (2sin(πω))^{2} with decreasing ω. Thus, the result concerning the low-frequency parts of spectra obtained from Eq.(3) applies also to that obtained from Eq.(1). For this reason, ${\scriptscriptstyle \frac{d(\cdot )}{dt}}$ is not further considered.
Fluctuations of a component x of an equilibrium solution at near-zero frequencies, that contribute to a non-zero Γ(0)-variance, cannot be determined by the fluctuations of its differential forcing f, since f de facto does not fluctuate at the lowest frequencies. But what else, if not f, generates the fluctuations of x at near-zero frequencies? Here we do not pretend to answer this question rigorously. Instead, we discuss a few ideas which we think are crucial for addressing this question. The overall outcome of the discussion – that the fluctuations of x at near-zero frequencies result from an integral effect – seems to be reasonable and self-evident, as these fluctuations definitely exist when the dynamical system of our concern is integrated forward in time.
We start by considering for a given dynamical system, the major frequencies (or the major timescales) of the underlying dynamics, defined as the frequencies (timescales), at which the spectra of f have their maxima. For the Lorenz model (Figures 7b and 8a), the most pronounced spectral peak is found for the differential forcing of the third component (blue line) located at the frequency near 1.42 × 10^{–2} cpts, which corresponds to a timescale of about 70 time steps. For the first and second component (green and red), the spectra of f are less peaky. There is a maximum located at frequencies slightly lower than 1.42 × 10^{–2} cpts and a maximum at frequencies around 2 × 10^{–2} cpts. We hence consider the frequencies of a few 10^{–2} cpts as the major frequencies of the Lorenz dynamics. The maxima of the spectra of x (Figure 7a) are located at similar frequencies as the spectra of f, although the spectra either further increase with decreasing frequency or drops only slightly at frequencies lower than the major frequencies. Consistent with these spectrum maxima, the time series of x and f (Figure 2) fluctuate visually only on the major timescales of the underlying dynamics. For x, whose spectrum has a white low-frequency plateau, the total low-frequency variance is evenly distributed over all low frequencies, making the fluctuations at an individual low frequency extremely weak and impossible to spot (even when sufficiently long records are used to resolve these low frequencies). Generally, we assume that the major timescale of a dynamical system of our concern is finite.
Given the finiteness of the major timescale of the underlying dynamics, it becomes even more striking that x fluctuates at zero and almost-zero frequencies. We approach this issue by considering a Brownian particle embedded in fluid molecules. The differential forcing of the Brownian particle describes the collisions of the particle with its surrounding fluid molecules. The major timescale of the underlying dynamics is the collision timescale. With each individual collision, the Brownian particle receives an impulse signal, which is a function of time that is only non-zero within a very short time interval and zero everywhere outside the interval. An impulse signal has an essentially white spectrum. Thus, a collision received by a Brownian particle makes the particle to fluctuate a little at almost all frequencies, especially the lowest frequencies.
The Brownian motion example can be extended to a dynamical system such as a GCM or the Lorenz model. A particular interaction of a component x with the other components, as described by the differential forcing f of x, produces a signal in x that is non-zero within a time interval, but zero outside the time interval. Since the length of the interval, which corresponds to the major timescale of the underlying dynamics, is always negligibly short relative to the infinite timescale, the signal can be considered as an impulse signal. Similar to a collision received by the Brownian particle, an interaction with other components received by x makes x to fluctuate a little at almost all frequencies, especially the lowest frequencies.
We emphasize that different from a delta function, which has a peak of zero duration, an impulse-like signal received by x has a “peak” lasting over a non-zero duration. Even for a collision of a Brownian particle with a molecule, which is completed within a time interval that is too short to be recognized by the human eye, the duration of the collision is, relative to the major timescale of the differential forcing that describes the collisions, not zero. Figure 9 shows typical spectra of impulse-like signals with a peak of non-zero duration. Different from the spectrum of a delta function, which is white and non-zero over the entire frequency domain, the spectrum of a signal of non-zero duration has a white low-frequency plateau. We emphasize also that x can only receive a full impulse-like signal, after x has completed the interaction with the other components, or in other words, after x has evolved along a solution over the duration of the interaction. Consequently, the low-frequency fluctuations arising from an impulse-like interaction received by x is not compatible with the low-frequency fluctuations of the differential forcing f, which is a function of every time instance during the interaction.
Along an equilibrium solution, component x receives not only just one interaction with other components of x, but many interactions. While each impulse-like interaction received by x makes x to fluctuate a little at the lowest frequencies, many impulse-like interactions can produce the non-zero Γ(0)-variance of x. We refer the net effect of many interactions of x with other components of x along a solution as the integral effect of all other components on component x. The net integral effect cannot be represented as a differential forcing defined at individual time instances, but can only be obtained by integrating f forward in time.
Still unclear is, as an equilibrium solution continues into an infinite time, what stops the integral effect to include more impulse-like interactions, thereby keeping the Γ(0)-variance and the respective total equilibrium variance finite. Purely integrating (i.e. simply adding) interactions of x with other components along a solution resembles a random walk that integrates a white noise over time. Such a pure integration makes the variance to increase with increasing time, and making an equilibrium state unsustainable. Obviously, there must be a dissipation, or a negative feedback as suggested by Hasselmann (1976), that counteracts the pure integration of interactions of x with other components.
Hasselmann identified the required dissipation as a part of the underlying dynamics, i. e. as a part of the differential forcing f. Such a dissipation dissipates variations of x at every time instance. We know from the example of Brownian motion that there is another type of dissipation: As a Brownian particle moves among fluid molecules, it experiences a friction that results from numerous collisions. This friction does not exist in the governing equations formulated in terms of differential forcing, since within a time interval δt → 0 that defines the differential operator ${\scriptscriptstyle \frac{d}{dt}}(\cdot )$ , a Brownian particle receives only isolating collisions. It appears that for a system in its equilibrium state reached under a constant external forcing, there exist two types of dissipation: One appears as a part of differential forcing and counteracts the constant external forcing at every time instance, which is needed to prevent the system to run away under the constant external forcing. The other appears as a part of the integral effect that does not exist in δt → 0 and counteracts the continuous integration of impulse-like interactions over time. As the impulse-like interactions make x to fluctuate at low frequencies, the latter type of dissipation is linked to (low-frequency) fluctuations, in line with the dissipation-fluctuation theorem (DFT) (Risken, 1989, chapter 7). DFT states that the response of a system to a change in its external forcing, which depends sensitively on the dissipation processes inside the system, is closely related to the natural fluctuations occurring without the change in the external forcing. Clearly, this dissipation issue and its link to DFT need further investigations.
Finally, we point out that the differential forcing f considered in this paper differs from a mechanical force f_{mech} considered in classical mechanics in two aspects. First, when exerting a mechanical force f_{mech} on an object, the force is equally received by every part of the object. For a system considered in this paper, although all components of x are subjected to the same external forcing, the differential forcings of any two different components are generally different. The external forcing is not equally received by all components of x. Secondly, a mechanical force f_{mech} exerted on an object can be long-lasting, whereas a differential forcing considered in this paper is always short living. In the simplest case in which f_{mech} is constant, the momentum of the object would be constantly accelerated, and the variance of the momentum would be infinite. In terms of spectral relation given in Eq.(33), which holds both for solutions of a system considered in this paper and for the evolution of the momentum caused by a mechanical force f_{mech}, the spectrum ${\Gamma}^{{f}_{mech}}(\omega )$ of f_{mech} would be a line spectrum with a non-zero value at ω = 0, and the spectrum of the momentum of the object would be proportional to ${\Gamma}^{{f}_{mech}}(\omega ){(2\pi \omega )}^{-2}$ and hence infinite at ω = 0. Different from that, a differential forcing considered in this paper does not have any energy (variance) at frequency zero, and de facto does not have any energy at all near-zero frequencies (see Figure 8a). Such a forcing is ephemeral. The ephemerality of f is reflected in decaying correlation functions, a feature not existing in classical mechanics.
The ephemerality is linked to the multidimensionality of the dynamical system considered in this paper. For a multidimensional system, even though the form of f(x) as a function of x does not change over time, the contributions of different components of x to f(x) generally change over time. In case of a Brownian motion, the contributions to f(x) at different times stem from collisions of the Brownian particle with different molecules. Because of this constant change in contributors to f(x), f cannot be long-lasting.
An ephemeral differential forcing f of x generates fluctuations in two different ways: one via f at each time instance; the other via the ephemerality of f related to the decaying and reawakening of f over time – due to the constant external forcing and the dissipation contained in f. The former makes x to fluctuates at the frequency of f, in the same way as in classical mechanics that the momentum variations at a frequency result from the variations of the mechanical force at the same frequency. The latter, which can only be realized by integrating the dynamical system forward over time, makes x to fluctuate at the lowest frequencies via impulse-like interactions of x with other components of x.
This paper considers a dynamical system described by a multidimensional state vector x such that the temporal evolution of a component x of x is governed by ${\scriptscriptstyle \frac{dx}{dt}}=f(\text{x}).$ . Under a constant external forcing, the system reaches an equilibrium state characterized by equilibrium fluctuations that have well-defined finite equilibrium variances and decaying covariance functions. We show that an equilibrium variance can only be decomposed in frequency domain when taking all equilibrium solutions of the considered equilibrium state into account. We show further that for a system of our concern, the differential operator ${\scriptscriptstyle \frac{d(\cdot )}{dt}}$ , together with the white low-frequency spectral plateau of x, constrain the differential forcing f of x such that f de facto does not fluctuate at near-zero frequencies. The variance related to the fluctuations of x at near-zero frequencies can hence not be generated by fluctuations of f. Instead, it results from an integral effect of many impulse-like interactions of x with other components of x occurring along an equilibrium solution over time. This integral effect cannot be represented as a differential forcing and only emerges when integrating the system forward in time. The forward integration sets the arrow of time toward an equilibrium state.
The significance of the present study depends on the amount of the total equilibrium variance that results from the integral effect. We assume that this variance equals the Γ(0)-variance, i.e. the variance that is associated with the white low-frequency spectral plateau of x. For the Lorenz system, we estimate the Γ(0)-variance in terms of ν defined as
ω_{*} is the frequency that for ω > ω_{*} the spectrum of component x, Γ^{x}(ω), deviates from its white low-frequency plateau. For the three Lorenz components, ω_{*} is indicated by the colored vertical lines in Figure 7a). ν amounts to 13.86%, 11.48%, 3.79% of the total variance, respectively, suggesting that only a small portion of the equilibrium variance is determined by the integral effect. We speculate that given a component x of a dynamical system, the amount of Γ(0)-variance depends on the “inertia” of the component. A more inert component needs to interact with other components a greater number of times before it starts to fluctuate. The integral effect that results from a larger number of interactions of x with other components leads to a higher white low-frequency plateau of x, and with that a larger Γ(0)-variance of x. It is conceivable that the Γ(0)-variance of a Brownian particle is essentially identical to its entire equilibrium variance.
This paper has strong implication for inert components of a dynamical system. The equilibrium variance of an inert component, which is largely determined by the integral effect, can no longer be deciphered by considering processes described by the differential forcing at individual time instances. If we are not content with estimating the total equilibrium variance by integrating the governing equations in form of ${\scriptscriptstyle \frac{dx}{dt}}=f(\text{x}).$ forward in time, but wish to cast the integral effect in an evolution equation, we need to find a new framework that quantifies and represents the integral effect. Such a framework can take the form of a stochastic difference equation, which represents the role of the integral effect in generating variance in terms of a stochastic forcing, and at the same time, contains a dissipation counteracting the continuous generation of variance by the stochastic forcing. Here I deliberate use the word “difference”, as oppose to “differential”, to emphasize that the integral effect is an effect not expressible in term of the differential operator ${\scriptscriptstyle \frac{d}{dt}}(\cdot )$ . Within an infinitesimal time interval associated with the operator ${\scriptscriptstyle \frac{d}{dt}}(\cdot )$ , the integral effect does not exist; we have only the differential forcing at every time instance. The integral effect is the physical “thing” brought up in Section 1 that escapes the underlying deterministic dynamics represented by the differential forcing f. In this sense, a framework that quantifies the integral effect based on stochastic difference equations is not an approximation of the underlying deterministic differential equation, but represents processes not represented by the differential forcing. The subject of such a framework are fluctuations at near-zero frequencies only, rather than fluctuations at the major frequencies of the underlying dynamics.
The concept of stochastic climate models proposed by Hasselmann (1976) is an attempt toward formalizing such a framework. If such a framework is established for the dynamical system of our concern (e.g. for the climate), we would have for any one component x of the state vector x two types of evolution equations: One describes the detailed interaction of x with other components occurring at every time instance in terms of the differential forcing f(x) of x. The other describes the generation and the dissipation of the portion of the equilibrium variance of x that arises from the integral effect occurring along a solution over time. The total equilibrium variance of x can be obtained by integrating the first type of evolution equations, but cannot be fully explained by the dynamics represented by only one of the two types of evolution equations.
All data and subsequently the figures are generated by Matlab-scripts archived in the publication repository of the Max Planck Society under http://hdl.handle.net/21.11116/0000-0009-9904-6.
The additional file for this article can be found as follows:
AppendicesAppendix A to D. DOI: https://doi.org/10.16993/tellusa.25.s1
I thank Christian Reick, Cathy Hohenegger, Eduardo Zorita, Hans von Storch and Jochem Marotzke for helpful discussions. My thanks go also to the two anonymous reviewers for their constructive comments.
The authors have no competing interests to declare.
Balaji, V, Maisonnave, E, Zadeh, N, Lawrence, BN, Biercamp, J, Fladrich, U, Wright, G, et al. 2017. Cpmip: Measurements of real computational performance of earth system models in cmip6. Geoscientific Model Development, 10(1): 19–34. Retrieved from https://gmd.copernicus.org/articles/10/19/2017/. DOI: https://doi.org/10.5194/gmd-10-19-2017
Bartlett, M. 1946. On the theoretical specification of sampling properties of auto correlated time series. J.R. Stat. Soc., B8: 27–41. DOI: https://doi.org/10.2307/2983611
Bloomfield, P. 1976. Fourier Analysis of Time Series: An Introduction. Wiley.
Dijkstra, HA and Ghil, M. 2005. Low-frequency variability of the large-scale ocean circulation: A dynamical systems approach. Reviews of Geophysics, 43: RG3002. DOI: https://doi.org/10.1029/2002RG000122
Eckmann, J and Ruelle, D. 1985. Ergodic theory of chaos and strange attractors. Reviews of Modern Physics, 57(3): 617–656. DOI: https://doi.org/10.1103/RevModPhys.57.617
Gardiner, C. 2004. Handbook of Stochastic Methods for physics, chemistry, and natural sciences. Springer. DOI: https://doi.org/10.1007/978-3-662-05389-8
Hasselmann, K. 1976. Stochastic climate models Part I. Theory. Tellus, 28(6): 473–485. DOI: https://doi.org/10.1111/j.2153-3490.1976.tb00696.x
James, S and James, PM. 1989. Ultra-low-frequency variability in a simple atmo spheric circulation model. Nature, 342: 53–55. DOI: https://doi.org/10.1038/342053a0
Jenkins, GM and Watts, DG. 1968. Spectral analysis and its applications. Holden-Day.
Kerr, RA. 2000. A North Atlantic climate pacemaker for the centuries. Science, 288: 1984–1985. DOI: https://doi.org/10.1126/science.288.5473.1984
Lasota, A and Mackey, M. 1994. Chaos, fractals, and noise: Stochastic aspects of dynamics. In Applied mathematical sciences, 97: 472. 2nd ed. New York: Springer. DOI: https://doi.org/10.1007/978-1-4612-4286-4
Lorenz, E. 1963. Deterministic nonperiodic flow. Journal of Atmospheric Sciences, 20(2): 130–141. DOI: https://doi.org/10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2
Mitchell, J. 1976. An overview of climatic variability and its causal mechanisms. Quaternary Research, 6(4): 481–493. Retrieved from https://www.sciencedirect.com/science/article/abs/pii/0033589476900211. DOI: https://doi.org/10.1016/0033-5894(76)90021-1
Newman, M, Alexander, MA, Ault, TR, Cobb, KM, Deser, C, Lorenzo, ED, Smith, CA, et al. 2016. The pacific decadal oscillation, revisited. Journal of Climate, 29(12): 4399–4427. Retrieved from https://journals.ametsoc.org/view/journals/clim/29/12/jcli-d-15-0508.1.xml. DOI: https://doi.org/10.1175/JCLI-D-15-0508.1
Panagiotopoulos, F, Shahgedanova, M and Stephenson, D. 2002. A review of Northern Hemisphere winter-time teleconnection patterns. Journal de Physics IV, 12: 27–47. DOI: https://doi.org/10.1051/jp4:20020450
Priestley, M. 1981. Spectral analysis and time series. Academic Press.
Rasmusson, E and Carpenter, T. 1982. Variations in tropical sea surface temperature and surface wind fields associated with the Southern Oscillation- El Niño. Monthly Weather Review, 110: 354–384. DOI: https://doi.org/10.1175/1520-0493(1982)110<0354:VITSST>2.0.CO;2
Risken, H. 1989. The Fokker-Planck equation. In (p. 472). Springer. DOI: https://doi.org/10.1007/978-3-642-61544-3
von Storch, H and Zwiers, FW. 1999. Statistical analysis in climate research. Cambridge University Press. DOI: https://doi.org/10.1007/978-3-662-03744-7_2
von Storch, J-S. 1999a. The reddest atmospheric modes and the forcings of the spectra of these Modes. Journal of the Atmospheric Sciences, 1614–1626. DOI: https://doi.org/10.1175/1520-0469(1999)056<1614:TRAMAT>2.0.CO;2
von Storch, J-S. 1999b. What determines the spectrum of a climate variable at zero frequency? Journal of Climate, 2124–2127. DOI: https://doi.org/10.1175/1520-0442(1999)012<2124:WDTSOA>2.0.CO;2