## 1. Introduction

Linear stability theory for fluid systems has been extensively studied because of its role in advancing understanding of physical phenomena, including structure and growth of perturbations, growth of errors in forecast models, transition from laminar to turbulent flow, and maintenance of the turbulent state. Historically, linear stability theory has addressed problems of deterministic growth using the method of modes (Rayleigh 1880; Charney 1947; Eady 1949). However, the method of modes is incomplete for understanding perturbation growth even for autonomous systems because the nonnormality of the linear operator in physical problems often produces transient development of a subset of perturbations that dominates the physically relevant growth processes (Kelvin 1887; Farrell 1982, hereafter F82). Recognition of the role of nonnormality in linear stability led to the development of generalized stability theory (GST; F82; Farrell 1988; Farrell and Ioannou 1996a, hereafter FI96a). Compared to the methods of modes, the methods of GST, which are based on the nonnormality of the linear operator, allow a far wider class of stability problems to be addressed, including perturbation growth associated with aperiodic time-dependent certain operators to which the method of modes does not apply (Farrell and Ioannou 1996b, hereafter FI96b; Farrell and Ioannou 1999, hereafter FI99). An example of such an aperiodic time-dependent stability problem is the forecast error growth problem in which the nonnormality of a certain but time-dependent linear system, the tangent linear operator of the forecast, produces asymptotic Lyapunov instability (FI99). In addition to addressing problems of deterministic growth of sure initial perturbations, GST also addresses problems of the growth of statistical distributions of perturbations to certain operators allowing prediction of statistical quantities including variance, fluxes, and structures in turbulent flow (Farrell and Ioannou 1993a,b), which has led to a new mechanistic model of geophysical and laboratory turbulence (Farrell and Ioannou 1994, 1995, 1998).

There remains a class of linear stability problems still to be addressed by the methods of GST, which is the stability of uncertain systems. The problems described above involve growth of perturbations in a system with no time dependence or a system with known time dependence: the perturbations to this system may be sure or stochastically distributed and may be imposed at the initial time, or distributed continuously in time, but the operator to which the forcing is applied is considered to be certain. However, it may happen that we do not have complete knowledge of the system that is being perturbed. For instance, the parameterizations of damping and radiation in the tangent linear forecast equations may not be certain but may rather include realistic statistical variability about their mean values. Another example is departures from the climatological mean of realizations of the planetary-scale Northern Hemisphere flow. These realizations are drawn from an ensemble of statistically similar flows but, as is well known, there are significant temporal variations of the planetary waves about the mean. In these examples the system governing the growth and structure of perturbations depends not only on the mean state but also on uncertain but statistically known variations of the state about its mean. Three regimes in the statistical analysis of uncertain system are distinguished: systems in which the time dependence of the statistical fluctuations of the operator are temporally correlated for short intervals compared to the damping and oscillation timescales of the associated mean operator, systems in which the statistical fluctuations of the mean operator are correlated for long time intervals compared to the timescales of the mean operator, and the physically important transitional case of operator fluctuation on timescales comparable to those of the mean operator.

Describing the stability of uncertain systems requires introducing the distinct concepts of sample stability, moment stability, and ensemble mean stability (Arnold 1984; Arnold and Kliemann 1987). Remarkably, uncertain systems that are sample (Lyapunov) stable may be higher-order moment (e.g., variance) unstable and uncertain systems that are sample and moment unstable may be ensemble mean stable. We show that a single realization of a time-dependent system exhibits only the sample (Lyapunov) asymptotic exponential growth for all its moments.

Compared to growth of ensemble mean quantities themselves, growth of ensemble mean second-moment quantities corresponds more closely to notions of stability and transient growth as generally understood. The sample (Lyapunov) growth rate is attained asymptotically by any single realization of the second moment, yet ensemble second-moment quantities can be unstable even for sample stable systems. We interpret these results physically using analysis based on the ensemble probability density function.

Ensemble mean and moment quantities are physically interpretable as applying to single realizations when the associated operator is stable and forced so that the ergodic assumption linking the ensemble and time mean can be made, in which case the ensemble mean and moment quantities apply to observations made using a time-averaging instrument. In Part I of this paper we show that in the ensemble mean of an uncertain system there is an optimal structure producing the greatest response to harmonic forcing and so the ensemble mean frequency response as well as the structure of the optimal forcing and response in the ensemble mean at each frequency apply to time averages of a sufficiently long time series from a single realization, if the ergodic assumption is valid.

In Farrell and Ioannou (2002, hereafter Part II) growth and stability of second-moment quantities are examined in more detail and the optimal perturbation problem for initial and stochastic forcing is cast and solved in the context of the second moment. We illustrate methods for analyzing ensemble mean growth in uncertain systems using simple examples including barotropic and baroclinic flow models with fluctuating winds and fluctuating damping parameters.

## 2. Moment stability in uncertain flows

*ψ*(

*t*) produced by realizations of an uncertain system is defined as

*ψ*(

*t*), produced by statistically similar realizations of a dynamical system. It is also useful (cf. Arnold 1984; Arnold and Kliemann 1987) to consider the asymptotic ensemble mean exponent of the first moment of the state vector,

*λ*

_{0}

*λ*

_{1}

*λ*

_{2}

*λ*

_{n}

*g*, takes the value of 2 or ½, with equal probability over unit time intervals. The Lyapunov exponent is

*p*-moment exponent, as defined, is given by

^{1}

We can use probability density functions of the ensemble member growth rates to examine moment growth more closely. If the growth rate of a perturbation is calculated after *n* units of time there are *n* + 1 possible growth rates distributed on the interval [log(1/2), log(2)] = [–0.69, 0.69] according to the binomial distribution (Fig. 1). As *n* increases the distribution of growth rates approaches a *δ* function about the Lyapunov exponent, *λ*_{0} = 0. However, despite the concentration of the distribution around *λ* = 0, for every *n* the growth rates remain distributed over the interval [–0.69, 0.69], and despite the vanishing probability as *n* → ∞ of any given realization achieving a growth rate significantly different from *λ*_{0}, the higher-moment exponents remain sufficiently influenced by the finite extent of this binomial distribution to give for all *n* the nonzero growth rates obtained in (9) (right panel of Fig. 1).

*P*

_{t}(

*λ*), of growth rates over an interval

*t*,

*p*moment exponents are

*λ*

_{0}≤

*λ*

_{1}≤

*λ*

_{2}≤ · · · ≤

*λ*

_{N}is clear from (13). In contrast to the finite higher-moment exponents of the binomial pdf (Fig. 1), Gaussian pdf's inevitably lead to instability of sufficiently high moments because the Gaussian can produce ensemble members, however rare, with unbounded growth rates that dominate the higher-moment exponent obtained in (13) (Arnold and Kliemann 1987). The pdf's of physical quantities are often approximately Gaussian near their mean but with vanishing probability of extreme values so that a distribution similar to the binomial may be more appropriate than the Gaussian in predicting higher-moment growth for physical systems. Moreover, for practical purposes the distinction between the Lyapunov exponent and moment exponents fades as the integration interval increases because it becomes increasingly unlikely that any realization deviates from the Lyapunov exponent and structure.

*t*→ ∞ the pdf becomes an infinitely narrow distribution about the Lyapunov exponent is known as the multiplicative ergodic theorem of Oseledec (Arnold 1998); this allows calculation of the Lyapunov exponent of an uncertain system from a single sufficiently long integration of the perturbation equations. On the other hand, as we have seen, the growth rate of higher-order moments depends on the finite extent of the growth rate pdf even as this pdf becomes increasingly centered about the Lyapunov exponent. If, for example, the distribution of growth rates

*P*

_{t}(

*λ*) is a Gaussian of width

*δ*(

*t*) = {〈

*λ*

^{2}(

*t*)〉}

^{1/2}, then according to (13) the Lyapunov exponent of the

*p*moment over the time interval

*t*is

*e*

^{pλt}in (14) is obtained by a cumulant expansion using the fact that the higher-order cumulants of a Gaussian distribution vanish (Van Kampen 1992):

*p*-moment exponent is then given exactly by

*t*→ ∞ asymptotic dependence of

*λ*

_{p}(∞) on the width of the pdf,

*δ*(

*t*), is clear. If we assume that

*δ*

^{2}(

*t*) =

*O*(1/

*t*), as when the growth rate distribution is Gaussian, then asymptotically the

*p*-moment Lyapunov exponent increases linearly in the moment

*p*, confirming that regardless of the value of the Lyapunov exponent,

*λ*

_{0}= lim

_{t→∞}〈

*λ*(

*t*)〉, sufficiently high-order moments will be unstable.

In order to illustrate these ideas consider the example of the modified Eady problem in which the zonal flow varies about a constant mean shear. The model is Boussinesq with constant stratification on an *f* plane and has periodic boundary conditions in the zonal, *x*, direction; solid walls in the meridional, *y*, direction; and a solid lid at height *z* = *H*. The zonal flow *U*(*z*, *t*, *ω*) is generalized to be time dependent and uncertain, a realization of which, *ω*, is taken from a probability space Ω. Horizontal scales are nondimensionalized by *L* = 1200 km, vertical scales by *H* = *fL*/*N* = 10 km, velocity by *U*_{0} = 50 m s^{–1}, and time by *T* = *L*/*U*_{0}, so that a time unit is approximately 6.7 h. The Brunt–Väisälä frequency is *N* = 10^{–2} s^{–1}, and the Coriolis parameter is *f* = 10^{–4} s^{–1}.

*ψ*for each realization

*ω*is

*ψ*(

*z*,

*t*,

*ω*)

*e*

^{ikx+ily}, where

*k*is the zonal and

*l*is the meridional wavenumber: the perturbation potential vorticity is

*D*

^{2}

*ψ*, with

*D*

^{2}≡ ∂

^{2}/∂

*z*

^{2}–

*k*

^{2}–

*l*

^{2}and the perturbation potential vorticity damping rate is

*r*.

*z*= 0, and tropopause,

*z*= 1, provides the boundary conditions

*U*′(0,

*t*,

*ω*) and

*U*′(1,

*t*,

*ω*) denote the velocity shear at

*z*= 0 and

*z*= 1, respectively. The coefficient of Ekman damping Γ

_{g}≡

*N*/

*U*

_{0}(

*ν*/2

*f*)

^{1/2}has value Γ

_{g}= 0.0632, corresponding to a vertical eddy momentum diffusion coefficient

*ν*= 20 m

^{2}s

^{–1}in the boundary layer.

*U*(

*z,*

*t,*

*ω*) =

*z*+

*u*(

*z,*

*t,*

*ω*) is composed of an ensemble mean component,

*z*, and uncertain deviations,

*u*(

*z*,

*t*,

*ω*). The linear perturbation dynamics is

**A**is the sure Eady operator,

**A**

*D*

^{2}

^{–1}

*ikzD*

^{2}

*rD*

^{2}

**A**

_{1}(

*t*,

*ω*) is the uncertain operator associated with the wind fluctuations,

*D*

^{2}has been rendered invertible by incorporating the boundary conditions.

*u*(

*z*,

*t*,

*ω*) =

*ϵ*Σ

_{i}

*ξ*

_{i}(

*t*,

*ω*)

*u*

_{i}(

*z*), in which the wind fluctuations result from modulation of a set of wind profiles,

*u*

_{i}(

*z*), weighted by the identically distributed random variable,

*ξ*

_{i}(

*t*,

*ω*), and

*ϵ*is an amplitude factor. In terms of the time-independent operators

**B**

_{i}the uncertain operator has the form

**B**

_{i}is the certain operator,

*η*

_{i}(

*t*,

*ω*) are generated by the stochastic differential equation,

*d*

*η*

_{i}

*ν*

*η*

_{i}

*dt*

*dW*

*W*a Wiener process. In (26) the symbol

*η*has been used in place of

*ξ*as an indication that the random variable is Gaussian. The standard deviation of

*η*

_{i}is (2

*ν*)

^{–1/2}and the autocorrelation time is 1/

*ν*. The wind profiles were chosen to give variance increasing with height, and also to almost always produce westerly winds. Typical realizations of the resulting mean flow with (〈

*ξ*

^{2}

_{i}

^{1/2}= 0.3 and autocorrelation time

*t*

_{c}= 1 are shown in Fig. 2.

Consider growth rates over the relatively short interval 8.5 days with meridional wavenumber *l* = 2, noise autocorrelation time 7 h, and an rms fluctuation amplitude *ϵ* = 0.3. Moment growth rates as a function of wavenumber are shown in Fig. 3. The middle curve is the ensemble mean 8.5-day growth rate as a function of wavenumber with the bars indicating the standard deviation from the ensemble mean. Individual realizations result in widely different growth rates over this time interval. The top curve shows the exponent of the second moment, *λ*_{2}, and demonstrates the ordering *λ*_{0} ≤ *λ*_{2}. For reference the optimal growth rate for the mean flow over the same interval is also shown. Because of the time dependence of the flow, the ensemble mean growth rate has increased for almost all wavenumbers.

The pdf of the growth rate obtained from 10^{5} simulations over *t* = 80 days is shown in Fig. 4. This time interval was chosen to be long enough so that initial conditions do not influence the result, and the mean of the pdf is a good approximation of the asymptotic Lyapunov exponent. In this example the Lyapunov exponent is negative but the second and higher moments are unstable, as shown in Fig. 5. The large *p* growth rates shown are not the asymptotic values because the limited number of simulations in this Monte Carlo calculation do not explore the unlikely events responsible for producing the higher-moment growth rates. If the number of simulations is increased, the higher-order moments asymptote to higher values, and in the limit of an infinite number of simulations the moment exponents grow linearly with *p*, as required by (16) given Gaussian wind fluctuations. However, as previously mentioned, the Gaussian pdf overestimates the probability of rare events in many physical problems so this asymptote is not likely to be observed.

Consider an ensemble of flows that are Lyapunov (sample) stable while at the same time being second-moment unstable. Such an ensemble would have divergent energy and variance although the ensemble mean magnitude of the perturbations decays to zero. This is because relatively rare high growth trajectories dominate the second moment (e.g., positive growth rates in the pdf shown in Fig. 4). Associated with these high growth rate realizations are structures that differ from the structures associated with the sample mean, which grows at the rate of the top Lyapunov exponent and has the structure of the first Lyapunov vector. The most unstable second-moment structure is the top EOF, which is calculated from the covariance matrix; EOFs for uncertain systems are examined in Part II of this paper.

## 3. Dynamics of the ensemble mean field

Consider a statistical description of perturbation growth in which the perturbation dynamics is governed by an uncertain linear operator with realizations 𝗔(*t*, *ω*), where *ω* is taken from a probability space Ω. We concentrate on uncertain operators 𝗔(*t*, *ω*) of the form 𝗔(*t*, *ω*) = 𝗔 + *ϵ**η*(*t*, *ω*)𝗕, with *η*(*t*, *ω*) a Gaussian process^{2} with zero mean, unit variance, and specified autocorrelation time *t*_{c}.

We have seen that the stability properties of an uncertain system can be understood from the pdf of the growth rate. For example, because the width of the growth rate pdf becomes vanishingly small with time the Lyapunov exponent and related structure can be obtained almost surely from analysis of a single long integration of the perturbation equations. Therefore, analysis of the structure of the first Lyapunov vector does not require consideration of ensembles (cf. FI99). We turn now to the problem of analyzing ensemble average growth and growth of higher-order moments. If we start with a sure initial condition in an (additively) unforced but multiplicatively uncertain system, a difference among the higher-order moment growth rates and the Lyapunov growth rate exists only in an ensemble of flow realizations and it seems at first that our study requires Monte Carlo simulation. However, this computationally expensive task can in many cases be avoided, because deterministic equations exist governing the evolution of higher-order moments in uncertain systems under appropriate restrictions.

An equation for the evolution of the ensemble average field, 〈*ψ*(*t*)〉, and the second moment, 〈*ψ*_{i}(*t*)*ψ*^{*}_{j}*t*)〉, can be readily obtained if the temporal variability of the fluctuations is taken to be white in time (Arnold 1992). Although it is a great advantage for analysis to take the temporal randomness in the uncertain operators to be white, this assumption is often hard to justify because white noise processes in common with Gaussian processes have unbounded fluctuations, for example, implying infinite wind speeds or negative damping rates. We wish to obtain equations that are not restricted in validity to the white noise limit, at least in part to assess the validity of the white noise limit.

We begin with derivation of an equation for the ensemble mean. While ensemble mean stability is not directly related to moment stability, it has a distinct and interesting interpretation that we will discuss, and it underpins the development of equations for the second moment. Uncertain second-moment dynamics, which is the subject of Part II of this paper, is of particular interest because it is associated with perturbation energy and the EOFs of the stochastically forced uncertain system, and its stability properties correspond more clearly to usual notions of perturbation growth.

We first derive the ensemble mean equation valid for the case in which the time mean operator 𝗔 commutes with the perturbation operator 𝗕, so that the commutator [𝗔, 𝗕] = 𝗔𝗕 – 𝗕𝗔 vanishes. We then deal with the case in which the commutator [𝗔, 𝗕] is nonzero.

*Dynamics of the ensemble mean field for commuting mean and fluctuation operators*

*t*) = 𝗔 +

*ϵ*

*η*(

*t*)𝗕, and in which the operators 𝗔 and 𝗕 are assumed to commute. Because they commute, the operators 𝗔 and 𝗕 can be simultaneously diagonalized by a change of coordinates and the evolution of the perturbations can be followed in these normal coordinates.

^{3}The behavior of this

*n*-dimensional system in normal coordinates is that of

*n*uncoupled one-dimensional uncertain systems. The exact equation for the mean can be obtained if

*η*(

*t*) is a Gaussian process with zero mean and unit variance that decorrelates with time at rate

*ν*, that is, 〈

*η*(

*t*

_{1})

*η*(

*t*

_{2})〉 = exp(–

*ν*‖

*t*

_{1}–

*t*

_{2}‖ ). For such a process the autocorrelation time,

*t*

_{c}, is

*t*

_{c}= 1/

*ν*. Because [𝗔, 𝗕] = 0, the propagator from

*t*to

*s*is

*θ*(

*t*,

*s*) ≡

^{t}

_{s}

*η*(

*s*)

*ds*. If the pdf of

*θ*,

*p*[

*θ*(

*t*,

*s*)], is known then the mean of the propagators is

*p*[

*θ*(

*t*,

*s*)], can be obtained for all times because the central limit theorem guarantees that the integral,

*θ*(

*t*,

*s*), of the zero mean Gaussian processes,

*η*(

*t*), will also be distributed Gaussian with pdf:

**I**is the identity. It can be easily verified using (32) that the mean of the propagators, 〈

**P**(

*t*,

*s*)〉, does not satisfy the semigroup property 〈

**P**(

*t*, 0)〉 = 〈

**P**(

*t*,

*s*)〉〈

**P**(

*s*, 0)〉, implying that the mean field at time

*t*resulting from a sure initial condition at time

*s*, 〈

*ψ*(

*t*)〉 = 〈

**P**(

*t*,

*s*)〉

*ψ*(

*s*), satisfies a differential equation depending explicitly on the initial time

*s*. In the particular case of evolution from

*t*= 0 the mean field 〈

*ψ*〉 ≡ 〈

*ψ*(

*t*, 0)〉 satisfies

*t*= 0 on this useful equation, which is exact in commuting case with Gaussian noise.

*ν*→ ∞) for all

*t*≫ 1/

*ν*we obtain asymptotically the ensemble mean equation

From (33) it is clear that for all autocorrelation times and all noise levels the asymptotic mean behavior of the uncertain system follows (34). Because this is the same ensemble mean behavior that would be obtained if the noise were white with variance 2*ϵ*^{2}/*ν*, the operator on the rhs of (34) will be referred to as the equivalent white noise operator. Equation (34) shows that asymptotically the effect of the fluctuations on the ensemble mean field is to modify the mean propagator by a term proportional to 𝗕^{2}, which for imaginary eigenvalues of 𝗕 contributes an additional damping to that of the mean operator, 𝗔. This extra damping results from phase averaging over the realizations of the ensemble. In contrast, real eigenvalues of 𝗕 increase the growth over that produced by the mean operator. This increased asymptotic growth results from the convexity of the exponential function in a manner analyzed in the context of time-dependent stability by FI99.

While (33) is exact for Gaussian distributed *θ*(*t*), the restriction on validity of (33) to Gaussian pdf's limits its utility for fluctuations with long autocorrelation times, *t*_{c} = 1/*ν*. We interpret this restriction as follows: each realization of the propagator (28) (with *s* = 0) contains a factor *e*^{ϵθ(t)B} in which *θ*(*t*) ≡ ^{t}_{0}*η*(*s*) *ds* for *η*(*t*) Gaussian with standard deviation 1. Consider replacing the Gaussian distribution *η* by a restricted distribution of fluctuations *ξ*(*t*) that is distributed Gaussian up to fluctuation magnitude Ξ_{0}, with zero probability of obtaining values ‖ *ξ*(*t*) ‖ > Ξ_{0}. Then for each such realization the term in the exponent, *ϵ**θ*(*t*), will be at most *ϵ*Ξ_{0}*t*; and when the averaging is performed the mean exponent must be smaller than *ϵ*Ξ_{0}*t* so that if restricted pdf's are to be modeled as Gaussian distributions it must be required according to (32) that *ϵ*Ξ_{0}*t* > *ϵ*^{2}*t*/*ν* for *t* ≫ *t*_{c}, which in turn requires Ξ_{0} > *ϵ*/*ν* = *ϵ**t*_{c} = *K*. The nondimensional number *K* = *ϵ**t*_{c} is called the Kubo number. We conclude that (33) is accurate for *ξ*(*t*) with restricted pdf only for Kubo numbers that do not exceed the maximum magnitude of the fluctuation, here taken to be Ξ_{0}. The Kubo number can be interpreted heuristically as the range in standard deviations over which rms fluctuations need to be distributed Gaussian for (33) to hold. Physical quantities such as wind or damping fluctuations seldom retain Gaussian statistics much beyond two standard deviations, so *K* = 2 is a practical upper bound for application of (33) to physical problems.

*ϵ*𝗕 and a minus state for which the propagator is 𝗔 –

*ϵ*𝗕. If the probability of residence in each state decays exponentially with time as

*e*

^{–t/tc}

*ψ*(

*t*)〉

_{±}, which denote the ensemble mean on the condition that at time

*t*the fluctuation is in the ± state, obey the equations

*ψ*(

*t*)〉 = 〈

*ψ*(

*t*)〉

_{+}+ 〈

*ψ*(

*t*)〉

_{–}is obtained. It is clear from (36a) and (36b) that in the limit of long correlation times, that is, for

*ν*= 1/

*t*

_{c}→ 0, the propagator of the ensemble mean becomes the average propagator of the two states without dependence on the correlation time. This prediction is in stark disagreement with the Gaussian ensemble mean propagator (35), which contains the mean fluctuation parameter

*ϵ*

^{2}/

*ν*that dominates the dynamics even for arbitrarily small

*ϵ*if fluctuations correlate for a sufficiently long interval (i.e.,

*ν*→ 0). We see that this behavior occurs because of the generally unphysical unbounded fluctuations permitted by a Gaussian pdf.

## 4. Ensemble mean dynamics for quasi-static fluctuations

*t*

_{c}= 1/

*ν*→ ∞ the rms amplitude of the operator fluctuation,

*ϵ*, tends to a finite nonzero value. In such cases the propagator can be approximated by its quasi-static limit in which the time dependence of the operator 𝗔 +

*ϵ*

*ξ*(

*t*)𝗕 is neglected but not its randomness. The ensemble mean propagator in the quasi-static approximation is consequently the propagator averaged over the fluctuation realizations

*P*(

*ξ*). The quasi-static approximation is formally valid for

*t*≪

*t*

_{c}. However, its validity extends for all times if all realizations of 𝗔 +

*ϵ*

*ξ*(

*t*) 𝗕 lead to perturbation decay, with decay times shorter than

*t*

_{c}. Note that (37) is also valid whether or not 𝗔 and 𝗕 commute and is not limited to Gaussian fluctuations.

*η*, which are Gaussian with zero mean and unit variance we replace

*ξ*by

*η*in (37) to obtain the mean propagator

*P*(

*ξ*):

*ϵ*

*ξ*𝗕 are stable.

## 5. A barotropic channel example in which the mean and fluctuation operator commute

*ψ*(

*t*)

*e*

^{ikx}, where

*x*is the coordinate in the zonal direction, are governed by the barotropic vorticity equation

*U*

_{0}is the mean flow,

*r*is the damping rate,

*ϵ*is the amplitude of the wind fluctuations, and the flow fluctuations

*η*(

*t*) are assumed to be Gaussian with autocorrelation time

*t*

_{c}= 1/

*ν*. In this case the mean operator 𝗔 = –

*ikU*

_{0}–

*r*and the fluctuation operator 𝗕 = –

*ik*commute and the ensemble mean equation (33) is

*t*≪

*t*

_{c}, (43) predicts Gaussian decay,

*t*≫

*t*

_{c}, (43) reduces to the white noise approximation,

*ψ*decays at the damping rate

*r*independent of wavenumber, (45) predicts that the ensemble mean, 〈

*ψ*〉, obtained by averaging over all realizations has an additional rate of decay that is asymptotically proportional to the square of the wavenumber, as if the ensemble were subject to a diffusion with coefficient

*ϵ*

^{2}/

*ν*. This effective diffusion results from cancellation among the realizations and is not a property of any single realization. This effective diffusion is a general feature of vacillating advection because the fluctuation operator 𝗕 is of the form

**u**· ∇, where

**u**is the fluctuating wind, and consequently the 𝗕

^{2}operator in (34) has the form of an anisotropic diffusion

**u**· ∇(

**u**· ∇) (cf. Knobloch 1978).

Consider a mean flow with *U*_{0} = 5, *r* = 0, and perturbations with wavenumber *k* = 1/5. The mean of 2000 realizations of (41) with these parameters is shown in Fig. 6. The solution found using the exact propagator (43) is indistinguishable from this mean. Remarkably, even for autocorrelation time *t*_{c} = 3 and fluctuation amplitude *ϵ* = 0.3 (*K* = 0.9) the asymptotic ensemble mean propagator (45) accurately captures the decay of the ensemble mean for all times (top panel). An example of the evolution of a sure initial condition for a larger Kubo number (*K* = 2) is shown in Fig. 6 (bottom panel); in this case the equivalent white noise prediction (45) overestimates the decay.

The Kubo number appropriate for wind fluctuations in this model is estimated from the expression *K* = *ku*′*t*_{c}, where *u*′ is the amplitude of the wind fluctuations, and *k* the perturbation wavenumber. A 30–m s^{–1} mean flow speed and a 1000-km perturbation spatial scale result in an advective timescale of 0.4 days. Flow fluctuations with 4-day autocorrelation and rms speed 1/10 that of the mean flow correspond to nondimensional autocorrelation time *t*_{c} = 10 and nondimensional fluctuation amplitude *ϵ* = 1/10, which gives Kubo number *K* = 1. Because we expect an approximate Gaussian pdf for wind fluctuations up to one standard deviation with *K* = 1, the Gaussian approximations is valid.

*t*≪

*t*

_{c}is

*t*≫

*t*

_{c}the limit is

## 6. Ergodic interpretation of the ensemble mean

The ensemble mean can be straightforwardly interpreted as the result of averaging over a large number of realizations as in taking the mean of an ensemble of forecasts. Alternatively, although the ensemble mean is not defined for a single realization, it is of theoretical and practical interest to determine under what circumstances the ensemble mean can be interpreted as applying to a single realization.

We make the ergodic assumption that if the system reaches a stationary state then its time mean is the same as its ensemble mean. This assumption has been verified using Monte Carlo simulations for uncertain systems excited by both additive deterministic and stochastic forcing. Note that identifying the ensemble mean with the time mean requires that the system be forced so that a nontrivial stationary state is obtained.

**I**is the identity. To arrive at the above expression we used (32) and assumed that the integrals exist, which requires that the operator 𝗔 +

*ϵ*

^{2}/

*ν*𝗕

^{2}be asymptotically stable. The ensemble average resolvent

**R**(

*ω*) that maps the sinusoidal forcing at frequency

*ω*to the ensemble average response at this frequency is then identified as

*ω*= 0. The above results require that 𝗔 and 𝗕 commute.

*ω*. If station measurements are made and Fourier analyzed to obtain the frequency response, then with fluctuations a weakening of the response and a broadening of the spectrum occurs, in accord with the theory presented above, with the observed resolvent in the case of a fluctuating wind being

The frequency response (52) is shown in Fig. 7 for fluctuation magnitude *ϵ* = 1/3 and two autocorrelation times, (left panel) *t*_{c} = 1/3 (*K* = 1/9) and (right panel) *t*_{c} = 6 (*K* = 2). For reference the response with no fluctuations (53), and the quasi-static response (40) are also shown. For short autocorrelation times the response is close to the mean response, while for *t*_{c} = 6 (*K* = 2) the response is almost identical to the quasi static. The exact frequency response in (52) has been derived on the assumption that the fluctuations are Gaussian, and it is not expected to describe in practice fluctuations which for large autocorrelation times have large corresponding Kubo numbers. As the Kubo number increases with increasing autocorrelation times the validity of (52) requires that the fluctuations retain their Gaussianity over a range of standard deviations equal to the Kubo number, *K*, something that is not expected in practice. For such cases the ensemble average frequency response is better approximated by the quasi-static approximation or by modeling the fluctuations by a telegraph Markov process as in appendix A.

In multidimensional systems the resolvent matrix contains information about the system's response to the spatial structure of the forcing. If the forcing structure is *F*, as in (49), then under the ergodic assumption the response at frequency, *ω*, as obtained by an averaging instrument, is **R**(*ω*)*F*, where **R**(*ω*) is the resolvent of the ensemble mean equation. The resolvent matrix contains full information about the response of the system to any forcing structure and at any frequency but this may be more information than we need and the maximum singular value of the resolvent and its associated right singular vector may be sufficient information. The right singular vector of **R**(*ω*) is the optimal forcing structure leading to the maximum response at frequency *ω* as recorded by an averaging instrument.

It should be noted that in uncertain systems the frequency response recorded by an averaging instrument is not the rms frequency response. The rms response, the square of which is the energy spectrum, is obtained from the resolvent of the second moment.

## 7. Ensemble mean dynamics for fluctuation operators that do not commute with the mean operator

*η*(

*t*) a Gaussian process with zero mean and unit variance that decorrelates exponentially in time at rate

*ν*.

*ϕ*(

*t*) by

*ψ*(

*t*) ≡

*e*

^{At}

*ϕ*(

*t*). The interaction perturbation obeys the equation

*ϕ*(0) =

*ψ*(0), the interaction perturbation at time

*t*is

*ϕ*

*t*

**G**

*t*

*ψ*

**G**(

*t*) is the time-ordered exponential

*η*(

*t*) this averaged propagator is exactly

*t*= 0 that is valid for all Kubo numbers:

*t*

_{c}≪

*t*

_{d}, where

*t*

_{d}= ‖ 𝗔 ‖

^{–1}is the decay time of the most persistent mode of the mean operator, we obtain from (63) the ensemble equation valid for general 𝗔 and 𝗕:

The ensemble mean equation (64) is also obtained in the limit of small Kubo number by perturbation techniques for non-Gaussian noise (Bourret 1962; Keller 1962; Papanicolaou and Keller 1971; Brissaud and Frisch 1974; Van Kampen 1974; see also Frisch 1968; Van Kampen 1992).

The white noise limit of (64) is obtained if the ratio *ϵ*^{2}/*ν* approaches a constant *D*, as *ϵ* → ∞ and *ν* → ∞. The rms of the white noise process is then *e*_{w} = (2*D*)^{1/2} and the augmented operator governing evolution of the ensemble mean is 𝗔 + *e*^{2}_{w}^{2}, which recovers the Stratonovich correction for white noise–induced drift in the Fokker–Planck equation (Arnold 1992).

Although computationally convenient, the white noise limit is not physically realistic, because it implies unbounded fluctuations. It is remarkable that the equivalent white noise approximation (64), with a different interpretation of *D* = *ϵ*^{2}/*ν*, is also valid for non–white noise processes with sufficiently short autocorrelation times. In examples even operators with a 2-day autocorrelation time are described well by (64).

If there is more than one independent structure so that an instantaneous realization of the total operator is 𝗔 + Σ_{i} *ϵ**ξ*_{i}(*t*)𝗕_{i}, with independent random *ξ*_{i}(*t*) having zero mean and unit variance, then the augmented deterministic operator that governs the evolution of the ensemble mean in (64) becomes 𝗔 + Σ_{i} *ϵ*^{2}_{i}*ν*^{2}_{i}

*Examples of ensemble mean dynamics for fluctuation operators that do not commute with the mean operator*

*x*and velocity

*υ*be governed by the equation

^{2}= 0, in the short autocorrelation time limit the mean position evolves as if there were no fluctuation in the spring constant [cf. (64)] and the ensemble mean equations are

This property of the dynamics of the harmonic oscillator may be better appreciated by examining Fig. 8, where it can be seen that the ensemble average of 500 realizations follows quite accurately the ensemble mean equation (66). This result is obtained only for short autocorrelation times for the spring constant fluctuations compared to the period of the oscillation. For such short autocorrelation times the position and the fluctuations are effectively uncorrelated, that is, 〈*x*(*t*)*ξ*(*t*)〉 = 0 as the modulation of the spring constant directly affects the acceleration, which only through integration affects the position of the oscillator. The effect of finite autocorrelation time is shown in Fig. 9. While for short autocorrelation times the ensemble mean displacement follows the ensemble mean equation (66), for finite autocorrelation times the ensemble mean fields decay as can be seen in the middle and bottom panel. It is clear that the stability of the ensemble mean fields is indicative neither of the stability of any single realization nor of the stability of higher moments, all of which are unstable in this example.

## 8. Conclusions

Developing methods for analyzing perturbation dynamics in uncertain flows is important for advancing understanding of both observations and forecast. Uncertainty may arise due to ensembling of actual realizations as in observations of Northern Hemisphere winter planetary waves or from a statistical description of a quantity as in a friction parameterization in a forecast model that accounts for variations about the mean dissipation value. In the first example interpreting an observation requires understanding ensemble statistics; in the second interpreting the effect of a statistically described parameter requires related methods.

The concept of moment stability is central to the analysis of uncertain systems and we have used the pdf of growth rates to interpret the relationship among sample and moment growths. Understanding the range of validity of analysis methods is also important and in this work we have identified the Kubo number as a practical measure for the required extent of Gaussianity in modeling uncertain dynamics.

A central concept in GST for uncertain systems is identifying the ensemble and time means through the ergodic assumption. This identification allows application of the powerful methods for ensemble analysis to time means of a single realization. Using this assumption a wide class of stationary statistics of a forced uncertain system can be obtained including the frequency response of the ensemble mean as recorded by an averaging instrument.

Optimal excitation plays a central role in GST and in Part I of this paper optimal forcing structures and related optimal response structures have been obtained as a function of forcing frequency for the ensemble mean response as recorded by an averaging instrument.

## Acknowledgments

The authors thank Ludwig Arnold, Prashant Sardeshmukh, and Cécile Penland for helpful discussions. This work was supported by NSF ATM-0123389 and by ONR N00014-99-1-0018.

## REFERENCES

Arnold, L., 1984: A formula connecting sample and moment stability of linear stochastic systems.

,*SIAM J. Appl. Math.***44****,**793–802.Arnold, L., 1992:

*Stochastic Differential Equations: Theory and Applications*. Krieger, 228 pp.Arnold, L., 1998:

*Random Dynamical Systems*. Springer-Verlag, 586 pp.Arnold, L., and W. Kliemann, 1987: Large deviations of linear stochastic differential equations.

*Stochastic Differential Systems*, H. J. Engelbert and W. Schmidt, Eds., Lecture Notes in Control and Information Sciences, Vol. 96, Springer-Verlag, 117–151.Bourret, R. C., 1962: Stochastically perturbed fields, with application to wave propagation in random media.

,*Nuovo Cimento***26****,**1–31.Brissaud, A., and U. Frisch, 1974: Solving linear stochastic differential equations.

,*J. Math. Phys.***15****,**524–534.Charney, J. G., 1947: On the dynamics of long waves in a baroclinic westerly current.

,*J. Meteor.***4****,**135–162.Eady, E. A., 1949: Long waves and cyclone waves.

,*Tellus***1****,**33–52.Farrell, B. F., 1982: The initial growth of disturbances in a baroclinic flow.

,*J. Atmos. Sci.***39****,**1663–1686.Farrell, B. F., 1988: Optimal excitation of neutral Rossby waves.

,*J. Atmos. Sci.***45****,**163–172.Farrell, B. F., and P. J. Ioannou, 1993a: Stochastic dynamics of baroclinic waves.

,*J. Atmos. Sci.***50****,**4044–4057.Farrell, B. F., 1993b: Stochastic forcing of the linearized Navier–Stokes equations.

,*Phys. Fluids***5A****,**2600–2609.Farrell, B. F., 1994: A theory for the statistical equilibrium energy and heat flux produced by transient baroclinic waves.

,*J. Atmos. Sci.***51****,**2685–2698.Farrell, B. F., 1995: Stochastic dynamics of the midlatitude atmospheric jet.

,*J. Atmos. Sci.***52****,**1642–1656.Farrell, B. F., 1996a: Generalized stability. Part I: Autonomous operators.

,*J. Atmos. Sci.***53****,**2025–2040.Farrell, B. F., 1996b: Generalized stability. Part II: Nonautonomous operators.

,*J. Atmos. Sci.***53****,**2041–2053.Farrell, B. F., 1998: Perturbation structure and spectra in turbulent channel flow.

,*Theor. Comput. Fluid Dyn.***11****,**215–227.Farrell, B. F., 1999: Perturbation growth and structure in time-dependent flows.

,*J. Atmos. Sci.***56****,**3622–3639.Farrell, B. F., 2002: Perturbation growth and structure in uncertain flows. Part II.

,*J. Atmos. Sci.***59****,**2647–2664.Frisch, U., 1968: Wave propagation in random media.

*Probabilistic Methods in Applied Mathematics,*Bharucha-Reid, Ed., Vol. 1, Academic Press, 75–198.Keller, J. B., 1962: Wave propagation in random media.

*Proc. 13th Symp. on Applied Mathematics*, Providence, RI, American Mathematics Society, 227–246.Kelvin, and Lord, 1887: On the stability of steady and of periodic fluid motions.

,*Philos. Mag.***23****,**459–464. 529–539.Knobloch, E., 1978: Turbulent diffusion of magnetic fields.

,*Astrophys. J.***225****,**1050–1057.Kubo, R., 1962: Generalized cumulant expansion method.

,*J. Phys. Soc. Japan***17****,**1100–1120.Palmer, T. N., 2001: A nonlinear dynamical perspective on model error: A proposal for non-local stochastic-dynamic parameterization in weather and climate prediction models.

,*Quart. J. Roy. Meteor. Soc.***127****,**279–304.Papanicolau, G., and J. B. Keller, 1971: Stochastic differential equations with applications to random harmonic oscillators and wave propagation in random media.

,*SIAM J. Appl. Math.***21****,**287–305.Rayleigh, and Lord, 1880: On the stability or instability of certain fluid motions.

,*Proc. London Math. Soc.***11****,**57–70.Sardeshmukh, P., C. Penland, and M. Newman, 2001: Rossby waves in a stochastically fluctuating medium.

*Progress in Probability,*P. Imkeller and J. S. von Storch, Eds., Vol. 49, Birkhauser Verlag, 369–384.Van Kampen, N. G., 1974: A cumulant expansion for stochastic differential equations.

,*Physica***74****,**239–247.Van Kampen, N. G., 1992:

*Stochastic Processes in Physics and Chemistry*. Elsevier, 465 pp.

## APPENDIX A

### Ensemble Mean Equation for Telegraph Fluctuations

*t*,

*ω*) = 𝗔 +

*ξ*(

*t*,

*ω*)𝗕, with

*ξ*(

*t*,

*ω*) taking only the values ±1 and that the jump from one value to the other follows a stationary Markov process known as a telegraph process (Van Kampen 1992), in which the conditional probability that

*ξ*has value

*ξ*

_{2}at time

*t*+ Δ

*t*given it had the value

*ξ*

_{1}at time

*t*depends only on the time interval Δ

*t*. It is physically reasonable to assume that this conditional probability,

*T*(

*ξ*

_{1}‖

*ξ*

_{2}; Δ

*t*), for sufficiently small Δ

*t*is

*T*

*ξ*

_{1}

*ξ*

_{2}

*t*

*δ*

_{ξ1,ξ2}

*w*

*ξ*

_{1}

*ξ*

_{2}

*t*

*w*(

*ξ*

_{1},

*ξ*

_{2}) from one state to the other is

*ξ*

*t*

*r*

*ξ*

*t*

*ν*

*τ*

*t*

_{c}= 1/

*ν*.

*ψ*(0) is determined from the ensemble average of the propagator

**P**(

*t*) by

*ψ*

*t*

**P**

*t*

*ψ*

**P**(

*t*) for a given realization of

*ξ*(

*t*) is

*t*

_{i}=

*i*

*δ*

*t*, and

*t*

_{n}=

*t*=

*n*

*δ*

*t*. The average propagator is the sum of all realizations of the propagator weighted by the probability

*p*that a specific realization (

*ξ*

_{1}, … ,

*ξ*

_{n}) of

*ξ*(

*t*) occurs, that is,

**P**(

*t*)〉

_{α}, the sum of the terms in (A9) for which

*ξ*

_{n}has the well-defined value

*α*at time

*t*. Then 〈

**P**(

*t*)〉 = 〈

**P**(

*t*)〉

_{+}+ 〈

**P**(

*t*)〉

_{–}. Because

*ξ*is Markov, 〈

**P**(

*t*+

*δ*

*t*)〉

_{±}can be expressed in terms of 〈

**P**(

*t*)〉

_{±}as

*δ*

*t*,

**I**+ (𝗔 ±

*ϵ*𝗕)

*δ*

*t*, and (A2), (A4) and (A5) used to obtain the differential equation governing the evolution of the mean propagator. Because the conditional average propagators 〈

**P**(

*t*)〉

_{±}determine the corresponding conditional average of the state 〈

*ψ*〉

_{±}, we immediately obtain the evolution equations for the state averages:

*ψ*〉 = 〈

*ψ*〉

_{+}+ 〈

*ψ*〉

_{–}obeys the exact mean equation,

This development is valid for noncommuting 𝗔 and 𝗕. The exact equation for the two state telegraphic process, (A14), or variants of it arise also under different noise assumptions, and (A14) can obtained directly from (A1) by perturbation techniques in the limit of small Kubo number, *K* ≪ 1 (Bourret 1962; Keller 1962; Brissaud and Frisch 1974; Van Kampen 1974).

## APPENDIX B

### Remarks on the Ensemble Mean Evolution Equation

*I*

_{n}=

^{t}

_{0}

*s*

^{n}

*e*

^{–νs}

*ds*. This equation is especially useful if the commutator vanishes at some order so that the series truncates at that order.

**G**(

*t*) may exponentially increase with time. This exponential growth of

**G**(

*t*) is bounded above by the difference between the decay rates of the most damped and the least damped mode of 𝗔. Because of the time dependence of

**G**(

*t*) and the fact that in general it does not commute with 𝗔, the exponential growth of

**G**(

*t*) does not necessarily imply the instability of the mean 〈

*ψ*〉. For example, consider the system

*η*a Gaussian noise process with zero mean, unit variance and autocorrelation time

*t*

_{c}= 1/

*ν*, and

**G**(

*t*) grows asymptotically as

*e*

^{(9–v)t}, but for

*ϵ*

^{2}/

*ν*= 1.33 and

*ν*= 1 Monte Carlo simulation confirms the prediction of (63) that 〈

*ψ*(

*t*)〉 decays as

*e*

^{–0.86t}. Note that the white noise limit in this case gives unstable growth at rate 0.33. Examples such as this in which

**G**diverges require exponentially small steps for accurate calculation of the propagator governing the ensemble mean evolution.

Realizations of the time-dependent velocity *U*(*z*, *t*) = *z* + *u*(*z*, *t*). The bold line is the mean wind, *U*(*z*) = *z*

Citation: Journal of the Atmospheric Sciences 59, 18; 10.1175/1520-0469(2002)059<2629:PGASIU>2.0.CO;2

Realizations of the time-dependent velocity *U*(*z*, *t*) = *z* + *u*(*z*, *t*). The bold line is the mean wind, *U*(*z*) = *z*

Citation: Journal of the Atmospheric Sciences 59, 18; 10.1175/1520-0469(2002)059<2629:PGASIU>2.0.CO;2

Realizations of the time-dependent velocity *U*(*z*, *t*) = *z* + *u*(*z*, *t*). The bold line is the mean wind, *U*(*z*) = *z*

Citation: Journal of the Atmospheric Sciences 59, 18; 10.1175/1520-0469(2002)059<2629:PGASIU>2.0.CO;2

Growth rates over 8.5 days as a function of wavenumber in the Eady model with uncertain fluctuating winds. The middle curve is the 8.5-day ensemble mean growth rate and the bars indicate one std dev from the ensemble mean value. The top curve shows the exponent of the second moment, which is by necessity larger than the mean growth rate. For reference the growth rate attained by the 8.5-day optimal perturbation on the mean flow is also shown. The time dependence of the flow leads to an increase of growth rate at almost all wavenumbers. The meridional wavenumber is *l* = 2, the coefficient of Ekman damping corresponds to a vertical diffusion coefficient of *ν* = 20 m^{2} s^{–1}, the autocorrelation time of the noise is 7 h, and the rms of the fluctuating components is 0.3

Growth rates over 8.5 days as a function of wavenumber in the Eady model with uncertain fluctuating winds. The middle curve is the 8.5-day ensemble mean growth rate and the bars indicate one std dev from the ensemble mean value. The top curve shows the exponent of the second moment, which is by necessity larger than the mean growth rate. For reference the growth rate attained by the 8.5-day optimal perturbation on the mean flow is also shown. The time dependence of the flow leads to an increase of growth rate at almost all wavenumbers. The meridional wavenumber is *l* = 2, the coefficient of Ekman damping corresponds to a vertical diffusion coefficient of *ν* = 20 m^{2} s^{–1}, the autocorrelation time of the noise is 7 h, and the rms of the fluctuating components is 0.3

Growth rates over 8.5 days as a function of wavenumber in the Eady model with uncertain fluctuating winds. The middle curve is the 8.5-day ensemble mean growth rate and the bars indicate one std dev from the ensemble mean value. The top curve shows the exponent of the second moment, which is by necessity larger than the mean growth rate. For reference the growth rate attained by the 8.5-day optimal perturbation on the mean flow is also shown. The time dependence of the flow leads to an increase of growth rate at almost all wavenumbers. The meridional wavenumber is *l* = 2, the coefficient of Ekman damping corresponds to a vertical diffusion coefficient of *ν* = 20 m^{2} s^{–1}, the autocorrelation time of the noise is 7 h, and the rms of the fluctuating components is 0.3

Pdf of the growth rate of perturbations over 80 days in the Eady model with fluctuating winds. The autocorrelation time of the fluctuations is *t*_{c} = 2 and their rms amplitude is *ϵ* = 0.6. The Lyapunov exponent is 〈*λ*〉 = –0.055; nevertheless, this flow is second-moment unstable with second-moment exponent 〈*λ*_{2}〉 = 0.05

Pdf of the growth rate of perturbations over 80 days in the Eady model with fluctuating winds. The autocorrelation time of the fluctuations is *t*_{c} = 2 and their rms amplitude is *ϵ* = 0.6. The Lyapunov exponent is 〈*λ*〉 = –0.055; nevertheless, this flow is second-moment unstable with second-moment exponent 〈*λ*_{2}〉 = 0.05

Pdf of the growth rate of perturbations over 80 days in the Eady model with fluctuating winds. The autocorrelation time of the fluctuations is *t*_{c} = 2 and their rms amplitude is *ϵ* = 0.6. The Lyapunov exponent is 〈*λ*〉 = –0.055; nevertheless, this flow is second-moment unstable with second-moment exponent 〈*λ*_{2}〉 = 0.05

Exponents of higher moments, *λ*_{p}, as a function of moment order, *p*. These exponents were obtained using the pdf shown in Fig. 4.

Exponents of higher moments, *λ*_{p}, as a function of moment order, *p*. These exponents were obtained using the pdf shown in Fig. 4.

Exponents of higher moments, *λ*_{p}, as a function of moment order, *p*. These exponents were obtained using the pdf shown in Fig. 4.

The real part of 〈*ψ*〉 obtained from a 2000-member ensemble of simulations of (41). Also plotted is the curve identical to that predicted from the exact propagator for the commuting case (43). The amplitude ‖ 〈*ψ*〉 ‖ according to (43) and the amplitude according to the equivalent white noise approximation (45) are indicated with the dotted line: (top) with autocorrelation time *t*_{c} = 3 and fluctuation amplitude *ϵ* = 0.3 (*K* = 0.9); the evolution follows closely the predictions of the white noise approximation; (bottom) with autocorrelation time *t*_{c} = 10 and fluctuation amplitude *ϵ* = 0.2 (*K* = 2); the equivalent white noise predictions (dot curve) overestimates the decay rate. The mean operator is 𝗔 = –*i* (the damping is *r* = 0) and the fluctuation operator is 𝗕 = –*i*/5.

The real part of 〈*ψ*〉 obtained from a 2000-member ensemble of simulations of (41). Also plotted is the curve identical to that predicted from the exact propagator for the commuting case (43). The amplitude ‖ 〈*ψ*〉 ‖ according to (43) and the amplitude according to the equivalent white noise approximation (45) are indicated with the dotted line: (top) with autocorrelation time *t*_{c} = 3 and fluctuation amplitude *ϵ* = 0.3 (*K* = 0.9); the evolution follows closely the predictions of the white noise approximation; (bottom) with autocorrelation time *t*_{c} = 10 and fluctuation amplitude *ϵ* = 0.2 (*K* = 2); the equivalent white noise predictions (dot curve) overestimates the decay rate. The mean operator is 𝗔 = –*i* (the damping is *r* = 0) and the fluctuation operator is 𝗕 = –*i*/5.

The real part of 〈*ψ*〉 obtained from a 2000-member ensemble of simulations of (41). Also plotted is the curve identical to that predicted from the exact propagator for the commuting case (43). The amplitude ‖ 〈*ψ*〉 ‖ according to (43) and the amplitude according to the equivalent white noise approximation (45) are indicated with the dotted line: (top) with autocorrelation time *t*_{c} = 3 and fluctuation amplitude *ϵ* = 0.3 (*K* = 0.9); the evolution follows closely the predictions of the white noise approximation; (bottom) with autocorrelation time *t*_{c} = 10 and fluctuation amplitude *ϵ* = 0.2 (*K* = 2); the equivalent white noise predictions (dot curve) overestimates the decay rate. The mean operator is 𝗔 = –*i* (the damping is *r* = 0) and the fluctuation operator is 𝗕 = –*i*/5.

Magnitude of the resolvent of the ensemble mean field as a function of frequency *ω*, in the barotropic flow example with fluctuating wind with rms amplitude *ϵ* = 1/3. Also shown are the quasi-static response (40) that gives the response for autocorrelation times that are very large, the response in the absence of fluctuations (53), and the exact response (52): (left) for autocorrelation time *t*_{c} = 1; (right) for autocorrelation time *t*_{c} = 6 in which case the exact response is very close to the quasi-static response. The mean wind is *U*_{0} = 5, the wavenumber is *k* = 1/5, and the coefficient of Ekman damping is *r* = 0.1. The fluctuation operator is 𝗕 = –*i*/5.

Magnitude of the resolvent of the ensemble mean field as a function of frequency *ω*, in the barotropic flow example with fluctuating wind with rms amplitude *ϵ* = 1/3. Also shown are the quasi-static response (40) that gives the response for autocorrelation times that are very large, the response in the absence of fluctuations (53), and the exact response (52): (left) for autocorrelation time *t*_{c} = 1; (right) for autocorrelation time *t*_{c} = 6 in which case the exact response is very close to the quasi-static response. The mean wind is *U*_{0} = 5, the wavenumber is *k* = 1/5, and the coefficient of Ekman damping is *r* = 0.1. The fluctuation operator is 𝗕 = –*i*/5.

Magnitude of the resolvent of the ensemble mean field as a function of frequency *ω*, in the barotropic flow example with fluctuating wind with rms amplitude *ϵ* = 1/3. Also shown are the quasi-static response (40) that gives the response for autocorrelation times that are very large, the response in the absence of fluctuations (53), and the exact response (52): (left) for autocorrelation time *t*_{c} = 1; (right) for autocorrelation time *t*_{c} = 6 in which case the exact response is very close to the quasi-static response. The mean wind is *U*_{0} = 5, the wavenumber is *k* = 1/5, and the coefficient of Ekman damping is *r* = 0.1. The fluctuation operator is 𝗕 = –*i*/5.

Ensemble mean displacement as a function of time for the harmonic oscillator with fluctuating spring constant. The initial condition is sure [*x*(0) = 1]. The ensemble mean displacement (thick line) oscillates as in the absence of fluctuations in the spring constant while individual realizations reveal the positive Lyapunov exponent resulting from the fluctuations. Two realizations are shown in the graph (thin lines). The noise is red with autocorrelation time *t*_{c} = 0.1 and rms amplitude *ϵ* = 0.6 giving Kubo number *K* = *t*_{c}*ϵ* = 0.06.

Ensemble mean displacement as a function of time for the harmonic oscillator with fluctuating spring constant. The initial condition is sure [*x*(0) = 1]. The ensemble mean displacement (thick line) oscillates as in the absence of fluctuations in the spring constant while individual realizations reveal the positive Lyapunov exponent resulting from the fluctuations. Two realizations are shown in the graph (thin lines). The noise is red with autocorrelation time *t*_{c} = 0.1 and rms amplitude *ϵ* = 0.6 giving Kubo number *K* = *t*_{c}*ϵ* = 0.06.

Ensemble mean displacement as a function of time for the harmonic oscillator with fluctuating spring constant. The initial condition is sure [*x*(0) = 1]. The ensemble mean displacement (thick line) oscillates as in the absence of fluctuations in the spring constant while individual realizations reveal the positive Lyapunov exponent resulting from the fluctuations. Two realizations are shown in the graph (thin lines). The noise is red with autocorrelation time *t*_{c} = 0.1 and rms amplitude *ϵ* = 0.6 giving Kubo number *K* = *t*_{c}*ϵ* = 0.06.

Ensemble mean displacement as a function of time for the harmonic oscillator with fluctuating spring constant. The rms amplitude of the fluctuation is *ϵ* = 1/7, and the natural frequency of the oscillator is *ω* = 1. (top) The ensemble mean displacement for *t*_{c} = 0.1; in this case there is negligible phase averaging and the solution is as predicted by the short autocorrelation time evolution equation (64), with the ensemble mean oscillating as if there were no fluctuations while the amplitude of individual realizations of the oscillator exponentially diverge. (middle) Intermediate autocorrelation time *t*_{c} = 1; the effects of correlation between displacement and velocity give a decay of the ensemble mean. (bottom) Long autocorrelation time *t*_{c} = 10 showing a large effect of phase averaging at this *t*_{c}.

Ensemble mean displacement as a function of time for the harmonic oscillator with fluctuating spring constant. The rms amplitude of the fluctuation is *ϵ* = 1/7, and the natural frequency of the oscillator is *ω* = 1. (top) The ensemble mean displacement for *t*_{c} = 0.1; in this case there is negligible phase averaging and the solution is as predicted by the short autocorrelation time evolution equation (64), with the ensemble mean oscillating as if there were no fluctuations while the amplitude of individual realizations of the oscillator exponentially diverge. (middle) Intermediate autocorrelation time *t*_{c} = 1; the effects of correlation between displacement and velocity give a decay of the ensemble mean. (bottom) Long autocorrelation time *t*_{c} = 10 showing a large effect of phase averaging at this *t*_{c}.

Ensemble mean displacement as a function of time for the harmonic oscillator with fluctuating spring constant. The rms amplitude of the fluctuation is *ϵ* = 1/7, and the natural frequency of the oscillator is *ω* = 1. (top) The ensemble mean displacement for *t*_{c} = 0.1; in this case there is negligible phase averaging and the solution is as predicted by the short autocorrelation time evolution equation (64), with the ensemble mean oscillating as if there were no fluctuations while the amplitude of individual realizations of the oscillator exponentially diverge. (middle) Intermediate autocorrelation time *t*_{c} = 1; the effects of correlation between displacement and velocity give a decay of the ensemble mean. (bottom) Long autocorrelation time *t*_{c} = 10 showing a large effect of phase averaging at this *t*_{c}.

^{1}

That the limit of *p* → 0 is the Lyapunov exponent *λ*_{0} has been shown in Arnold (1984) and Arnold and Kliemann (1987).

^{2}

The *ω* will henceforth be omitted for brevity.

^{3}

This use of the term normal refers to the orthogonality of the modes of 𝗔 in these coordinates and is unrelated to Gaussian normal distributions.