# Rough parameter dependence in climate models and the role of Ruelle-Pollicott resonances

See allHide authors and affiliations

Contributed by James C. McWilliams, November 22, 2013 (sent for review August 9, 2013)

## Significance

It is shown on a geophysical fluid model that the sensitive dependence on parametric variations of a given set of observations is related to the dominant modes of variability of the flow, as encoded through the Ruelle-Pollicott resonances. The latter are estimated by Markov modeling techniques. Such modeling and analysis methods can be used for any type of high-dimensional dissipative and chaotic dynamical systems.

## Abstract

Despite the importance of uncertainties encountered in climate model simulations, the fundamental mechanisms at the origin of sensitive behavior of long-term model statistics remain unclear. Variability of turbulent flows in the atmosphere and oceans exhibits recurrent large-scale patterns. These patterns, while evolving irregularly in time, manifest characteristic frequencies across a large range of time scales, from intraseasonal through interdecadal. Based on modern spectral theory of chaotic and dissipative dynamical systems, the associated low-frequency variability may be formulated in terms of Ruelle-Pollicott (RP) resonances. RP resonances encode information on the nonlinear dynamics of the system, and an approach for estimating them—as filtered through an observable of the system—is proposed. This approach relies on an appropriate Markov representation of the dynamics associated with a given observable. It is shown that, within this representation, the spectral gap—defined as the distance between the subdominant RP resonance and the unit circle—plays a major role in the roughness of parameter dependences. The model statistics are the most sensitive for the smallest spectral gaps; such small gaps turn out to correspond to regimes where the low-frequency variability is more pronounced, whereas autocorrelations decay more slowly. The present approach is applied to analyze the rough parameter dependence encountered in key statistics of an El-Niño–Southern Oscillation model of intermediate complexity. Theoretical arguments, however, strongly suggest that such links between model sensitivity and the decay of correlation properties are not limited to this particular model and could hold much more generally.

Sensitive behavior of long-term general circulation model (GCM) statistics is attracting increased attention (1⇓–3), but its origin and fundamental mechanisms remain unclear. These sensitive-behavior issues are of practical, as well as theoretical, importance in climate dynamics and elsewhere (4). For some GCMs, involving millions of variables, circumstances have been found where certain climate observables vary smoothly through a plausible parameter range (5) or where linear response theory applies over some range (6). On the other hand, this may not hold for every observable or parameter, and concerns arise regarding the role of some type of “structural instability” in sensitive parameter dependence (1, 2, 4).

The low-order Lorenz (L63) model (7) illustrates some of the relevant issues. Various statistics exhibit linear dependence over a broad range of parameters for which the dynamics is chaotic (e.g., figure 2 of ref. 8). The statistics’ linear dependence coexists here with structural instability of this model’s global attractor, as small variations in the parameters cause a plethora of topological changes (9). In particular, the unstable periodic orbits that appear or disappear as a parameter changes may only have a negligible effect on the model’s physical invariant measure (see below), if their period is longer than the decorrelation time of the dynamics.

In general, the role played by a system’s mixing and harmonic properties on the nature of its response has been only partially addressed. Only very specific results exist, in the deterministic setting, to support the idea that linear response of the long-term statistics (and of local variations of physical measures) may still hold in the absence of (topological) structural stability (10). For stochastic systems, more general results have been obtained (11), but it is still a challenge to relate the size of the parameter interval over which linear response may hold to the system’s mixing properties. Conditions for smooth but nonlinear response (e.g., quadratic) or else for rough parameter dependence—with many highly local variations in response over a given parameter interval—to occur are also poorly known.

To help us understand the circumstances in which one may expect one type of behavior rather than the other, we cast here this problem in a theoretical framework based on the modern spectral theory of dynamical systems (10, 12⇓⇓⇓⇓⇓⇓–19). The approach is illustrated on an El Niño–Southern Oscillation (ENSO) model of intermediate complexity. The model is governed by a system of coupled partial differential equations (PDEs), and it exhibits different degrees of roughness in its parameter dependence in different regimes. The relationship of statistics such as the power spectrum to dynamical features known as Ruelle-Pollicott (RP) resonances (20⇓–22) is outlined below and suggests the usefulness of estimating these resonances—despite the challenge of doing so in high-dimensional systems. To do so, we introduce here a unique approach that estimates these resonances as filtered through an observable chosen from the simulated scalar time series. This approach allows us to shed light on subtle relationships between the nonlinear mixing rate in the system’s phase flow and the nature of the parameter dependence of its long-term statistics.

## Intermediate ENSO Model and Its Key Properties

The intermediate-complexity ENSO model examined in this study is the Jin-Neelin (JN) model (23) forced by the seasonal cycle. The way we include the latter differs from the one used in ref. 24, but the main dynamical features are preserved (*SI Text*). The dynamics and thermodynamics of the resulting forced JN (fJN) model are based on those of the coupled ocean-atmosphere model of Cane and Zebiak (CZ) (25).

The fJN model’s main ingredients are the following. Its oceanic component is made up of two parts. The vertical-mean motions above the thermocline are governed by linearized shallow-water equations—forced by the wind stress—on an equatorial β-plane. The resulting currents drive an advection equation that describes the sea-surface temperature (SST) field at the Earth’s equator. The atmospheric component is a Gill-type model for the wind-stress anomaly field, which establishes a diagnostic relation (i.e., one with no time derivative present) between the latter and the SST anomalies. The magnitude of the wind stress anomalies controls the coupling between the oceanic and atmospheric components (*SI Text*).

We consider here, following ref. 23, a standard truncated version of this model, summarized in *SI Text*. The resulting numerical version has slightly more than 400 degrees of freedom: a modest number compared with a GCM but still a challenge for the study of statistical properties within the framework of transfer operator theory as presented elsewhere (12, 15, 17, 18, 26, 27).

The nonlinear interaction of the seasonal cycle in the fJN model with internal variability leads to a rich variety of dynamical behavior, from frequency-locked regimes to chaotic ones via a quasi-periodic route (24); the latter recalls the overlapping of Arnol’d tongues occurring in the CZ model (28). The internal variability arises through Hopf bifurcation when the basic state is steady (23), whereas in the fJN model, it arises via destabilization of a basic cycle of a 1-y period.

For a fixed coupling between the oceanic and atmospheric components, we analyze the response of various statistics to changes of a key model parameter *δ*. This parameter affects the travel time of the equatorially trapped waves (23, 24) that play an essential role in ENSO dynamics; it can yield significant changes in the Floquet spectrum of the linearized model (23). Here we examine regimes at strong nonlinearity over a modest range of *δ* (roughly 7% change) that nonetheless exhibit strong changes in various statistics of the simulated temperatures averaged over a Central Pacific region along the equator referred to as Niño-3.

## Parameter Dependence in the fJN Model

We begin by reporting two distinct types of parameter dependence (Figs. 1 and 2) with respect to small changes of the same primary parameter *δ*. For two different regimes, long model runs are produced to obtain these results over a fine *δ* grid: for each with and , an 8,800-y-long run is generated and sampled every 2 wk. The distinction between the two types of regimes is monitored by ; this parameter varies from zero to unity, and it controls the intensity of the anomalous surface-layer currents as a function of the wind stress anomalies (*SI Text*). When is close to unity, i.e., in the case of strong surface-layer feedback, stronger vertical and advection anomalies add to the rate of SST change (24).

Fig. 2 illustrates the *δ* dependence of the power spectrum of the Niño-3 SSTs. In the lower panel, for , a clear broad peak corresponds to a quasi-quadriennial (QQ) oscillation around 0.25 cycles/y (29) that occurs for most *δ* values. The upper panel, for , does not exhibit such a pronounced, broad interannual peak over the *δ* interval of study. For reasons that will become obvious later, we call the regime that corresponds to rapidly mixing and the one that corresponds to slowly mixing. For the moment, we can roughly say that these attributes are chosen in agreement with the decay rate of the autocorrelation function (ACF) of , which is typically faster for than for .

In both regimes and for , Figs. 1 and 2 illustrate the presence of *δ* intervals where chaos occurs, interleaved with intervals of periodicity, in agreement with the Arnol’d-tongue scenario noted above (24). Fig. 1 *B* and *D* demonstrates that in the slowly mixing case and on chaotic subintervals, sudden changes are manifested in the statistical moments of ; these changes may be relatively large within chaotic regimes (cyan dots in Fig. 1). At the same time, Fig. 1 *A* and *C* shows that nearly linear or, at least, smooth *δ* dependence takes place in the presence of chaos for the rapidly mixing case, except from jumps at transitions to periodic regimes. An exception occurs in a small range near the relative *δ* change value of 1 in the SD (Fig. 1*A*); this will be clarified in the penultimate section.

A more detailed inspection of the *δ* dependence of the statistics shown in Fig. 1 reveals interesting similarities with the *δ* dependence of the power spectrum plotted in Fig. 2. In the rapidly mixing case, the *δ* dependence of the power spectrum shown in Fig. 2*A* is rather smooth for the chaotic regions, as observed in Fig. 1 *A* and *C*, where the (narrow band) peaks are slightly modulated in magnitude across the *δ* interval of study. In the slowly mixing case—where the broad-band QQ peak is strongly modulated in magnitude and characterized by a rough *δ* dependence—striking correspondences are found between the changes of the SD shown in Fig. 1*B* and those of the QQ peak magnitude shown in Fig. 2*B*. This correspondence is natural given that the QQ peak captures much of the variance in this case.

Given this illustration of rough and not-so-rough parameter dependence within chaotic regimes, one is led to ask whether dynamical systems theory may help us clarify the relationship to broad-band energetic peaks and the rate of decay of correlations associated with these. Obviously this question is not limited to the fJN model and concerns chaotic dynamical systems in general. We recall first the relevant elements of the theory of RP resonances before presenting the main new ingredients, i.e., the Markov representations that will be used to provide concrete steps toward answering this question.

## RP Resonances and Decay of Correlations

In this section, we give a brief introduction to the spectral theory of dissipative dynamical systems (DDSs) (10, 12⇓⇓⇓⇓⇓⇓–19) while focusing on chaotic behavior. The next sections will show that this theory—when combined with appropriate Markov representations—provides a powerful set of concepts and tools that help us understand and quantify the occurrence of rough parameter dependence in complex systems, as illustrated in the previous section.

To simplify the presentation, let be a Euclidean vector space of dimension *d*, subject to a one-parameter group of smooth transformations , associated with the flow of a smooth but nonlinear system of ordinary differential equations, given by . The main objective of the spectral theory of chaotic DDSs is to study the evolution of probability laws induced by *S*_{t} instead of studying individual trajectories that exhibit chaotic behavior.

This objective is achieved by examining the family of Perron–Frobenius operators. Such a family , also known as transfer operators (12, 15, 17, 18, 26, 27), acts on probability measures ν and it is given byfor any measurable set . It gives the measure with respect to ν of the ensemble of points in *X* that occupy *E* at time *t*. Note that sometimes is denoted by , i.e., the pushforward of the measure *ν* by *S*_{t}.

Under mild assumptions on **F** and *ν*, it can be shown that is in fact a weak solution η emanating from ν at *t* = 0, in the Schwartz sense of distributions, of the transport equationon , where the operator div is the divergence operator on *X*. In parallel, the study of the evolution of densities associated with ν with due attention to the proper functional spaces in which these densities live (10, 12⇓⇓⇓⇓⇓⇓–19), is of prime importance for the theory.

It can be proven for hyperbolic dynamical systems and it is observed experimentally for many others (15, 30, 31) that a common feature of DDSs is the transformation of the initial Lebesgue measure into a measure that has still a finite density with respect to but that exhibits finer and finer structure, as time *t* evolves. Asymptotically, an invariant measure *μ* of Sinai-Ruelle-Bowen (SRB) type is generally reached (15, 30) as . This measure is physical in the sense that for , almost all , and any sufficiently smooth observable . This property is often referred to as the chaotic hypothesis that, roughly speaking, expresses an extension of the ergodic hypothesis to non-Hamiltonian systems (31).

Such a measure *μ* is typically singular with respect to while exhibiting smooth density in the expanding directions, or unstable manifolds, whose Haussdorf dimension is strictly less than *d*. In other words, the initial measure flows into microscopic scales of vanishing volume, whereas on the macroscopic scale, *μ*, supported by the unstable manifolds, is the only information from that remains visible after the dynamics acted over an infinite amount of time.

For Anosov flows (14), one considers functional spaces that capture such stretching and contracting effects of the dynamics and can then prove that the SRB measure *μ* may be equivalently characterized as the stationary solution of Eq. **2**. In ref. 14, it is then shown that the related statistical properties of the flow are accurately described by the spectral properties of the transfer operators acting on . Moreover, is a strongly continuous semigroup in , uniformly bounded in *t*, and its generator is given by , where **F** is the vector field generating the Anosov flow. As a consequence, the spectrum of is contained in the left-half complex plane, , and its resolvent —which determines the spectral properties of —is a well-defined bounded operator on that admits, for all and with , the following integral representation:

It can then be proven that the spectrum of the generator on consists only of isolated eigenvalues of finite multiplicity within a strip , for some that depends on the stretching and contracting rate of the dynamics (14); the rest of the spectrum is continuous and located in . In addition, the eigenspace associated with the null eigenvalue is spanned by the set of SRB measures; this set is reduced to only one such measure provided the zero is simple. We assume henceforth that the latter is the case, and *μ* will always refer to this unique SRB measure.

The extension of such rigorous results to other classes of chaotic DDSs is still a major challenge. However, a widespread conjecture is that the global picture is relevant to most chaotic DDSs. In other words, the spectrum of on some appropriate functional spaces for such a system consists always of the disjoint union of a continuous part and a discrete part. These two are called, respectively, the essential spectrum , and the point spectrum (13). When the existence of a (unique) SRB measure *μ* is assumed, classical function spaces, such as or , are suitable (17, 18).

In the interesting cases, i.e., when the point spectrum is not trivially reduced to , it is typically located in a vertical strip of the complex plane whose points are known as the RP resonances (20⇓–22). These resonances give precise information on correlation decay and on the power spectrum. We describe this information via formal mathematical arguments, referring to the specialized literature for their rigorous treatment.

Note first that, by making the change of variables where denotes the product map . This equation results from the general change-of-variables formula , with (32). Here the action of on the density *f* with respect to *μ* is defined by , the Radon-Nykodim derivative of with respect to *μ* (17, 18) (see also *SI Text*).

If we assume, without loss of generality, that *f* and *g* have a vanishing *μ*-ensemble mean, i.e., , then the right side of Eq. **4** denotes the correlation function . From the physical property of the SRB measure *μ*, this function equals for almost all , the more familiar cross-correlation coefficient at lag *t*, given by .

Using Eq. **3**, one can rewrite the Fourier transform of , given by , asThe meromorphic extension into the complex plane of , via Eq. **5**, tells us that the poles of the resolvent —which correspond to the RP resonances—introduce singularities into the complex Fourier transform, where is a complex frequency. These poles will manifest themselves in the power spectrum as peaks that stand out over a continuous background at frequency *ξ*, if the corresponding RP resonances with imaginary part *ξ* are close enough to the imaginary axis. The continuous background has its origin in the continuous part of .

Note that the position of the poles depends only on the system considered, whereas their residues depend on the observable monitored (21). As a result, some peaks in the power spectrum may disappear from observable to observable, depending on how large or small their residues may be. Ruelle (21, 22) and Pollicott (20) introduced this description of RP resonances as poles in the meromorphic extension of that are responsible for bumps in the power spectrum; they also connected these resonances to the decay rate of correlations.

The latter connections, between RP resonances and correlation decay, are subtler and closely related to the model flow’s relaxation toward the SRB measure *μ* as . If 0 is the only eigenvalue of on the imaginary axis, the system is mixing, and the existence of a gap between 0 and the rest of the spectrum governs the rate of convergence of toward *μ*. The size *τ* of this gap is given by .

More precisely, in such a case, one can prove that, for any probability measure *ν* that has density *ψ* with respect to the Lebesgue measure and for all suitably chosen test functions *φ*

The main step in establishing Eq. **6** is to prove the existence of the spectral gap ; remaining steps rely on the inversion of Eq. **3** and the properties of the resolvent 20. Furthermore, because and , one obtains from Eq. **6** that at the exponential rate *τ*. Thus, the spectral gap *τ* controls the decay of correlations, when the system is initiated out of the SRB equilibrium *μ*.

These arguments can be made rigorous in the context of Anosov flows (19), but, in general, the decay of correlations can be subexponential. This situation arises when the RP resonances are arbitrarily close to the imaginary axis. If , and the discrete spectrum is nontrivial, RP resonances lead to modulations in the decay of correlations, which correspond to peaks in the power spectrum.

Thus, RP resonances provide powerful theoretical tools to describe the variability in time of the system’s flow, in terms of spectral properties of the operator , where **F** is the nonlinear vector field that generates the flow. Key features of this variability include peaks in the power spectrum and the decay rate of correlations.

The original treatment of RP resonances (20, 22) was based, in fact, on a different approach, which used Markov partitions of the dynamics (15). The spectrum of —which lives within the unit disk—was analyzed, rather than that of , in the left-half plane. However, similar results relating the RP resonances to the decay of correlations and to the power spectrum were established. We preferred to follow here the framework of ref. 14 because of its connections with the fundamental Eq. **2**, whose analysis may benefit from PDE techniques (33) and the spectral theory of semigroups and their generators.

The latter theory offers conceptual advantages for describing the relations between the RP resonances and the decay of correlations on one hand and the power spectrum on the other. Nevertheless, for practical purposes, the approximation by a discrete time Markov process from a sequence of observations of the dynamics may be used to provide estimates of the (filtered) RP resonances when for instance the direct computation of is out of reach for large systems such as considered in this article. Such an approach is described below.

## Estimating RP Resonances from Observables

In a low-dimensional phase space, Markov partitions provide natural tools to study the spectral properties of or , the latter being approximated by a stochastic matrix *P* in the case of maps (26, 27). In essence, Ulam-type methods approximate the dynamics in phase space by a Markov chain whose transition probabilities are estimated from many simultaneous iterations of the map of interest over a large ensemble of initial data (26, 27). This approach can be rigorously justified for a large class of expanding or Anosov maps (27), and it can be used in the numerical estimation of RP resonances for low-dimensional models (34), although some drawbacks may arise in applications (35).

In a large-dimensional phase space, such as the one where the fJN model’s dynamics takes place, with , the Ulam approach becomes computationally intractable. A cheaper and less ambitious approach consists of taking a single observable *h* and, instead of trying to approximate the full transfer operator , seek a decomposition of the autocorrelation function associated with *h* from a long simulation into a sum of complex exponentials. In principle, this gives the positions of the RP resonances corresponding to nonvanishing residues associated with *h*. Padé approximation or Prony’s method are typically used (36). The main drawback of these techniques lies in the number of exponentials to be fitted to the signal, which may lead to an inaccurate estimation of RP resonance (35).

We propose here an intermediate approach based on sufficiently long model runs that exploit Ulam’s ideas but apply them to a Markov operator (17) that acts on a space of functions that depends only on the observed variables. As we will show, this operator is rigorously associated with the full transfer operator of the dynamics, given an observable *h* and the physical measure *μ* of the underlying map. The operator can then be approximated by a stochastic matrix *P*, which is estimated by computing a classical maximum likelihood estimator (MLE) (37), , from the sequence of observations .

When this sequence is long enough, the eigenvalues of provide crucial information about, and actual estimates of, the dominant RP resonances, as filtered by *h*, i.e., the resonances that correspond to nonvanishing residues and yield the largest contributions to the power spectrum of *h* (*SI Text*). Clearly, the more delicate point is the existence of such an operator , which we show hereafter for a broad class of chaotic systems. Details of a more practical nature about the approximation of then follow and are applied to the fJN model. Finally, note that our approach is complementary to the one in ref. 38 for Hamiltonian systems; the latter focuses on the related but distinct question of identifying metastable features of the dynamics.

### Markov Operators from Observables of Chaotic Systems.

In this section, we consider discrete dynamical systems given by , where , and **u**_{n} is assumed to be a periodic forcing of period *m*. We assume furthermore that —the solution operator that evolves the system from its state at time *p* to its state at time *n*—is well defined for any and that the semigroup property holds for all .

Chaotic behavior is synonymous here to the existence of a unique time-dependent, necessarily *m*-periodic SRB measure , which attracts the Lebesgue measure in a pullback sense, i.e., as , where *ρ* is the Lebegue measure (39). The latter property is assumed to hold in the absence of chaos as well, e.g., for a global limit cycle or for a quasi-periodic behavior. Note that in the presence of a positive Lyapunov exponent, due to a theorem of Ledrappier and Young, one can ensure the existence of such SRB measures by perturbing the governing evolution equations with an appropriate noise of very small intensity (39).

For simplicity, we denote now by *S* the map whose iterations give the states of the system at any integer multiple of the period *m* of the forcing. This map is a called the time-*m* map. By using an analog of Eq. **1** for discrete time, a transfer operator ℒ can be associated to *S*. By *m*-periodicity and from our assumptions on , the dynamical system generated now by *S* possesses as an SRB measure. We consider the system in this statistical equilibrium and define acting on densities *f* with respect to *μ* by for any such that .

Let be the (compact) support of *μ* and define the transition probability for the map *S* of reaching the Borel set *D* of from the Borel set *C* of bywhere is the characteristic function of set *A*, and the latter equality results from Eq. **4**, where for and . We now state the main result on which we will rely to explain the puzzling parameter dependence of the fJN model pointed out in Fig. 1 (*SI Text*).

### Theorem A.

*Let* *be a continuous observable of the dynamical system generated by* S, *with* . *Assume that* S *possesses a unique physical measure μ with support* . *Let* *be the set* *and* *be the pushforward of μ by h*. *Then there exists a Markov operator* *acting on* *such that* *and such that, for any Borel sets E and F of*

The proof of this result is a consequence of the general disintegration theorem of measures, from which can be constructed explicitly (*SI Text*). From Eq. **7**, we see that determines the transition probabilities for any pair (*C*, *D*) of Borel sets, and conversely, Eq. **8** shows that, when restricting the transition probabilities to pairs of the form , the Markov operator determines these exactly.

In general, the latter pairs run across a coarser family of subsets than the former, being a sub-*σ*-algebra of the Borel sets. We may thus say that characterizes a coarse-graining—induced by the observable *h*—of the actual dynamics, along with its transitions. At the same time, the operator given by theorem A offers a natural way to represent rigorously the sequence of observations as a finite-size sample of the discrete-time Markov process associated with (*SI Text*, corollary B). In brief, theorem A provides a rigorous basis for the assertion that the simple fact of observing a deterministic system allows us to represent the unobserved variables as noise (*SI Text*). The theory of Markov process can thus shed considerable light on the spectral properties of filtered by *h*, as we illustrate in the next sections.

### Stochastic Matrices from Observations and RP Resonances.

We assume hereafter that an SRB measure *μ* exists for the time *m* map *S* associated with the truncated version of the fJN model in with . The goal here is to analyze the spectrum of in and, more specifically, its dominant contributions. To approximate the dominant part of , we use a Galerkin procedure on a uniform grid of the one-dimensional set .

The eigenvalue problem is projected onto the problem in the subspace spanned by . The entries of the matrix are given by , where the second equality relies on theorem A. Thus, contains information about the actual transfer operator associated with *S*, as induced by *h*, and the supplementary coarse-graining induced by the partitioning of .

Note that, because is a Markov operator, its Galerkin approximation is a row stochastic matrix whose eigenvalues *λ* satisfy . In the Hamiltonian case of ref. 38, is self-adjoint and therefore—as *M* increases, i.e., the discretization becomes finer—one obtains immediately that the eigenvalues of approximate those of . For dissipative dynamics, however, is typically not self-adjoint, and only the dominant eigenvalues of can be robustly approximated (13, 16). Fortunately, it is the latter we are interested in here, and one can argue that the robustness of the dominant part of the spectrum does apply to the Markov operator associated with the fJN model (and *h*) (*SI Text*).

For a given *M*, a classical MLE, , can be used to approximate from a given sequence of observations , with being the Niño-3 SSTs. The entries of are then simply given by the relative frequencies , which converge to as with an error of order (37).

The dynamical interpretation of the Markov operator given by Eq. **8** and the numerical procedure described above provide a general framework for the estimation from a time series of the spectral gap in the RP resonances, as filtered by a particular observable *h*. This spectral RP gap can be estimated from the approximations of the dominant eigenvalues of , whenever the resolution *M* of the range of the observations *x*_{n} and the number *N* of observations are large enough; the statistical errors in these approximations can be estimated by bootstrapping arguments (*SI Text*). We describe next how the spectral RP gap turns out to be an important factor in the sensitivity of the empirical probability measure that approximates .

## Spectral RP Gap and Sensitivity of Statistical Equilibria

An important result in the stability theory of Markov chains was the discovery of sensitivity bounds relating the stability of a chain and its speed of convergence to equilibrium. Going back to the Markov representation of the dynamics of the fJN model given by theorem A, we take advantage of the sensitivity bounds for Markov chains to deduce sensitivity properties of the one-dimensional measure that the observable *h* extracts from the multidimensional SRB measure *μ*. Out of the many sensitivity bounds in the literature, it is the ones in ref. 40 that are most relevant here.

Recall first that the dual of defines a transition kernel (41), , which in turn defines a Markov operator on measures *ν* on given by , for any set *B* in the *σ-*algebra of the observed range of interest . The operators and are also linked by , with for (*SI Text*)*.* Uniform ergodicity will be the key concept here; it means that the iterates of the Dirac measures converge—uniformly in *x*—to in total variation (TV), i.e., there exists and such that, for all and all , .

We now consider perturbations of the Markov operator assumed to obey uniform ergodicity. The main result of ref. 40 stipulates that, if uniform ergodicity holds, then the invariant measure , associated with the Markov operator , satisfies the sensitivity boundwhere is the operator norm associated with the total variation norm on measures in (40). In Eq. **9**, is the smallest integer greater than or equal to , and (40).

The smallest *ρ* for which geometric ergodicity holds is called the rate of mixing of the Markov chain. For any fixed , grows superlinearly to infinity as , allowing in principle a large difference between and , even for perturbations that are relatively small, as measured by The dependence on *C* is much weaker, with increasing with *C*, but at a sublinear rate. The size of is thus strongly controlled by .

These results from ref. 40, together with theorem A, allow us to state—for the map *S* with SRB measure *μ*, as considered in the previous section—that if the associated Markov operator is uniformly ergodic, then the slower the mixing rate of , i.e., the closer is to 1, the larger we may expect the sensitivity of to be to perturbations of the system. From the dynamic interpretation of RP resonances, we conclude that regimes corresponding to slow decay with pronounced modulations—when observed through a given observable *h*—favor rough parameter dependence for the statistics built on *h*. These theoretical predictions are confirmed for the fJN model by the numerical calculations that follow.

## Spectral RP Gap for the fJN Model and Sensitivity

Recall that when the state space is finite and 1 is the unique, simple eigenvalue of a stochastic matrix *P* on the unit circle, then the mixing rate appearing in Eq. **9** is equal to , the subdominant eigenvalue of *P* (42). When *μ* is mixing, it can be shown that is irreducible and aperiodic, which in turn makes uniformly ergodic (27) such that Eq. **9** can be applied to . We adopt the strategy described above to estimate the gap between the unit circle and the subdominant eigenvalue of and to provide a confidence interval associated with this estimate for each value of *δ*.

Fig. 3 illustrates the use of the spectral gap to quantify the weak and strong mixing regimes discussed earlier. When the gaps observed in the chaotic regimes are compared between the case (*Lower*) and (*Upper*), we find— for each *δ* of interest—that the gap is typically smaller in the latter case than in the former. As a result, a higher sensitivity of the statistics is expected to occur in the case according to theorem A and Eq. **9**. Recalling the results of Fig. 1, we see that the numerical results in Fig. 3 along with their theoretical interpretations are in very good agreement with the experiments where the highest roughness in the *δ* dependence was observed for the case . Even the small regime of sudden changes in SD noted near relative δ-values of about 1 in Fig. 1*A* is consistent with the local decrease in the gap observed in Fig. 3 (*Bottom* panel), allowing higher sensitivity to occur locally in δ. Confidence intervals for these results are provided in the *SI Text*, supporting their robustness.

From the combination of the theory and numerical results, we infer that the occurrence of rough parameter dependence in the slow mixing case, where the QQ mode is the most energetic, is not a coincidence. This case corresponds to RP resonances that are closer to the imaginary axis than in the rapid mixing case. According to the sensitivity bound of Eq. **9** applied to the Markov representation of the dynamics provided by theorem A, this offers a favorable ground for sensitive behavior to occur.

## Concluding Remarks

This study is a first step in understanding the relationship between the time variability of a dissipative chaotic system and the parameter dependence of its long-term statistics. The relationship between these two aspects via the theory of RP resonances—and the data-based Markov representation developed in this article—opens up a wide range of possible investigations.

In this respect, other interesting Markov representations that have been used in climate dynamics (43⇓⇓⇓⇓⇓–49) might benefit from the framework of RP resonances, given their natural connection with the underlying nonlinear dynamics. In particular, applying RP resonances to the interpretation of metastability and flow regimes in connection to the LFV observed in geophysical flows (1, 44, 45, 50) and their possible role in parameter sensitivity, as well as in linear response theory in the presence of noise (11), are fascinating areas for further exploration.

## Acknowledgments

M.D.C. thanks Honghu Liu for help in preparing Fig. S1 and Joyce Meyerson for help in preparing Figs. 1–3. This study was supported by National Science Foundation Grants DMS 1049253 and AGS-1102838, Office of Naval Research Grant N00014-12-1-0911, and Department of Energy Grant DE-SC0006739.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: jcm{at}atmos.ucla.edu.

Author contributions: M.D.C., J.D.N., J.C.M., and M.G. designed research; M.D.C. and D.K. performed research; M.D.C. contributed new reagents/analytic tools; M.D.C. and D.K. analyzed data; and M.D.C., J.D.N., J.C.M., and M.G. wrote the paper.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1321816111/-/DCSupplemental.

Freely available online through the PNAS open access option.

## References

- ↵
- ↵
- McWilliams JC

- ↵
- Knight CG,
- et al.

- ↵
- ↵
- Neelin JD,
- Bracco A,
- Luo H,
- McWilliams JC,
- Meyerson JE

- ↵
- ↵
- ↵
- ↵
- Sparrow C

- ↵
- ↵
- ↵
- Baladi V

- ↵
- ↵
- ↵
- Collet P,
- Eckmann J-P

- ↵
- Keller G,
- Liverani C

- ↵
- Lasota A,
- Mackey MC

- ↵
- Sarig O

- ↵
- ↵
- ↵
- ↵
- ↵
- Jin F-F,
- Neelin JD

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Halmos PR

- ↵
- ↵
- ↵
- Horvat M,
- Veble G

- ↵
- ↵
- Billingsley P

- ↵
- ↵
- ↵
- ↵
- Kallenberg O

- ↵
- Hennion H,
- Hervé L

- ↵
- ↵
- ↵
- Ghil M,
- Robertson AW

- ↵
- ↵
- ↵
- ↵
- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Applied Mathematics

## Sign up for Article Alerts

## Jump to section

- Article
- Abstract
- Intermediate ENSO Model and Its Key Properties
- Parameter Dependence in the fJN Model
- RP Resonances and Decay of Correlations
- Estimating RP Resonances from Observables
- Spectral RP Gap and Sensitivity of Statistical Equilibria
- Spectral RP Gap for the fJN Model and Sensitivity
- Concluding Remarks
- Acknowledgments
- Footnotes
- References

- Figures & SI
- Info & Metrics