## Abstract

The combined effects of advection and diffusion on the equilibrium spatial structure of a tracer whose spatial variation is maintained by a large-scale forcing are considered. Motivated by the lower stratosphere, the flow is taken to be large-scale, time-dependent, and purely horizontal but varying in the vertical, with the vertical shear much larger than horizontal velocity gradients. As a result, the ratio *α* between horizontal and vertical tracer scales is large. (For the lower stratospheric surf zone *α* has been shown to be about 250.) The diffusion parameterizes the mixing effects of small-scale processes.

The three space dimensions and the large range between the forcing scale and the diffusive scale mean that direct numerical simulation would be prohibitively expensive for this problem. Instead, an ensemble approach is used that takes advantage of the separation between the large scale of the flow and the small scale of the tracer distribution. This approach, which has previously been used in theoretical investigations of two-dimensional flows, provides an efficient technique to derive statistical properties of the tracer distributions such as horizontal-wavenumber spectrum.

First, the authors consider random-strain models in which the velocity gradient experienced by a fluid parcel is modeled by a random process. The results show the expected *k*^{−1} Batchelor spectrum at large scales, with a deviation from this form at a scale that is larger by a factor *α* than the diffusive scale found in the absence of vertical shear. This effect may be crudely captured by replacing the diffusivity *κ* by an &ldquo=uivalent diffusivity” *α*^{2}*κ,* but the diffusive dissipation is then substantially overestimated, and the spectrum at large *k* is too steep. This may be attributed to the failure of the equivalent diffusivity to capture the variability of the vertical shear.

The technique is then applied to lower-stratospheric velocity fields. For realistic values of the diffusivity *κ,* the spectrum is found to be affected by diffusion at surprisingly large scales. For the value *κ* ∼ 10^{−2} m^{2} s^{−1} suggested in several recent papers, the diffusion is sufficiently strong that there is no clear *k*^{−1} regime, consistent with observations. The spectrum is then relatively well approximated by the observed *k*^{−2} power law in the range 20–200 km, but significantly steeper at smaller scales. For the molecular value *κ* = 10^{−4} m^{2} s^{−1}, in contrast, an unrealistic *k*^{−1} regime appears.

## 1. Introduction

The stretching properties of a flow are an important guide to its mixing ability. It is now recognized that effective stretching, and hence mixing, takes place not only in flows with a complex spatial and temporal structure, such as three-dimensional turbulent flows, but also in flows with simple spatial and temporal structure, which may be described as chaotic-advection flows. A key feature of chaotic-advection flows (and a useful working definition) is that stretching on small scales is dominated by the large-scale velocity. The corresponding regime of advection is called the “Batchelor regime” by the theoretical fluid-dynamics community. Such flows stand in strong contrast to three-dimensional turbulence, where it is the small-scale velocity field that dominates the small-scale stretching.

There is now much quantitative information on the stretching properties of stratospheric flows, obtained through studies based on velocity fields from large-scale meteorological datasets and calculating diagnostics such as contour lengthening rates or finite-time Liapunov exponents (e.g., Pierrehumbert 1991; Waugh et al. 1994; Chen 1994). These studies rely, implicitly or explicitly, on the assumption that chaotic advection (or the Batchelor regime) is relevant to stratospheric flows: if small-scale velocities dominated the stretching there would be little relevant information in the large-scale fields provided by meteorological datasets. The validity of the assumption has been tested by comparing results using winds at different spatial resolutions (Waugh and Plumb 1994; Methven and Hoskins 1999) and is also supported by basic fluid-dynamical arguments applying to strongly stratified, rapidly rotating flows (e.g., Bartello 2000; Shepherd et al. 2000).

However, quantification of stretching rates is not all that is needed to understand atmospheric chemical tracer distributions. For long-lived chemical tracers, the spatial structure of such distributions arises from the interplay between stretching on the one hand and molecular diffusion, perhaps enhanced by intermittent three-dimensional turbulence, on the other hand. Unless diffusion is neglected completely, the relation between stretching properties and spatial structure in tracer distributions is far from transparent, and it has recently been the focus of much interest in theoretical turbulence research (e.g., Shraiman and Siggia 2000; Falkovich et al. 2001).

One traditional diagnostic of the spatial structure of a chemical tracer (or of any other scalar quantity) is the power spectrum, *F*(*k*) say, regarded as a function of the modulus of the horizontal wavenumber *k* = |**k**|. Theoretical predictions are available for the form of the tracer spectrum in certain classes of flows. Comparing spectra derived from atmospheric observations with these theoretical predictions can therefore potentially give insight into atmospheric flows and atmospheric mixing processes.

In this paper, we focus on the tracer spectrum that is produced by large-scale quasi-two-dimensional winds and discuss theoretical predictions in relation to spectra derived from aircraft measurements of N_{2}O and O_{3} in the winter stratosphere (Strahan and Mahlman 1994; Bacmeister et al. 1996; Sparling and Bacmeister 2001). These measurements have revealed a clear power law for the spectrum in a wide range of scales (50–500 km, say), with an exponent close to −2. At smaller scales, the behavior is less certain: when raw data measured along the aircraft flight track are used, the power-law spectrum extends down to scales as small as 1 km (Bacmeister et al. 1996); however, this may be a result of vertical displacements of isentropic surfaces due, for example, to gravity waves rather than something that can be explained by large-scale advection. This was suggested by Sparling and Bacmeister (2001) who used a conditioning technique to remove the effect of vertical displacements. As a result, they obtain a spectrum that steepens significantly at small scales. The precise scale of transition appears to vary from year to year, and the small-scale spectrum is characterized by a power-law exponent of about −2.5.

The observed spectrum is difficult to reconcile with the idea that stratospheric flows are of the chaotic-advection type; that is, the velocity field, or more precisely the velocity gradients, are dominated by the large scales. (In terms of the kinetic energy spectrum this requires that the spectrum is *k*^{−3} or steeper.) Under this condition, and for a forced-dissipative situation in which a tracer is forced at large scales and diffused or dissipated at small scales, theory predicts a steady-state *k*^{−1} tracer spectrum—the Batchelor spectrum—in the so-called convective range of scales, between the forcing scale and the dissipation scale (e.g., Pierrehumbert 1994; Antonsen et al. 1996).

In an attempt to reconcile the spectra derived from measurements with the view of the stratosphere as dominated by large-scale winds, Ngan and Shepherd (1997) noted that a *k*^{−2} spectrum would be obtained if the tracer distribution included isolated edgelike structures associated with the edge of the polar vortex or with isolated filaments resulting from the occasional strong deformation of the vortex edge. They supported their arguments with contour-advection calculations, following the evolution of a piecewise-constant tracer field for a few days and showing that the resulting tracer spectrum varied as *k*^{−2}. In these simulations the isolation of edgelike structures and filaments was due to the finite integration time, though Ngan and Shepherd (1997) argued that in the real atmosphere filaments remain isolated because of mixing processes that limit their lifetime.

It seems to us that any complete explanation of observed stratospheric tracer spectra must be based on a model that includes both forcing and dissipation. For example, it is natural to regard the latitudinal tracer contrast on a given isentropic surface as continually forced on the large scale by the action of the meridional diabatic circulation, with the tracer distribution that is observed arising as an equilibrium between that large-scale forcing and the stirring effect of the flow. If one accepts this, then the situation studied by Ngan and Shepherd (1997) is not obviously relevant, and the process leading to a *k*^{−2} rather than a *k*^{−1} spectrum needs further explanation. Small-scale mixing processes may be a key part of the explanation since it is not clear on which range of scales they can be safely neglected. The purpose of the present paper is precisely to present an analysis of the effect of small-scale mixing processes, represented by diffusion, on the establishment of tracer spectra and, in particular, on tracer spectra in the stratosphere.

The computation of stationary tracer spectra for a large range of scales requires solution, with high spatial resolution and over a long time, of a three-dimensional forced advection–diffusion equation. This makes standard numerical techniques expensive to employ and calls for an alternative. In the chaotic-advection regime, because of the potential scale separation between tracer fields and velocity fields, such an alternative exists. Antonsen et al. (1996) use a technique based on a WKB approximation to solve the advection–diffusion equation, reducing the computation of the spectrum to the integration of ordinary differential equations for tracer wavenumber and amplitude along fluid-parcel trajectories. We employ this technique, which we describe in section 2. Antonsen et al. (1996) have successfully tested their technique in the case of two-dimensional randomized flows by comparing results with those obtained from explicit numerical simulations using a spectral code and, in doing so, verified the *k*^{−1} spectrum. Stratospheric flow can be regarded within a good approximation as two-dimensional but, as pointed out by Haynes and Anglade (1997), the layerwise nature of the flow and the associated vertical shear are crucial for tracer diffusion. The technique of Antonsen et al. (1996) is easily adapted to take that effect into account.

We study the impact of diffusion on tracer spectra using two types of models for the velocity field: theoretical random-strain models in section 3, and large-scale winds derived from observations in section 4. The ideas of chaotic advection suggest that random-strain models, which are based on the assumption that the strain experienced along a fluid-parcel trajectory is a random process, are relevant to stratospheric mixing (at least in the surf zone, where transport barriers have little influence). We do not attempt to confirm this in any precise way; we regard random-strain models as convenient tools to gain insight in the dependence of the spectrum on important parameters such as typical strain and shear, diffusivity, etc. This is particularly fruitful in two parameter limits for which simple analytical results are available. Our computations in section 4 based on observed winds are then interpreted in the light of the random-strain results.

Particular attention is paid to the effect of the vertical shear on the spectrum. In stretching and shearing flows such as the ones under consideration, the horizontal and vertical tracer scales decrease rapidly. On average, they decrease exponentially, at the same rate, and thus have a characteristic aspect ratio (Haynes and Anglade 1997). This suggests that the influence of the vertical shear can be modeled using an equivalent horizontal diffusivity: the shear is neglected, but the diffusivity *κ* is replaced by the larger value *κ*_{e} = *κ*(1 + *α*^{2}), where *α* is the characteristic ratio between horizontal and vertical scales (Haynes and Anglade 1997). Using both random-strain models and observed velocity fields, we examine the extent to which use of an equivalent diffusivity can lead to satisfactory predictions for the tracer spectrum. The shape of the spectrum in the diffusive range as well as the extent of this range depend strongly on the details of the vertical shear, in particular on its variability; a simple increase of the value of the diffusivity cannot therefore capture at all wavenumbers the effect of the vertical shear. In particular we find that, compared with the purely two-dimensional case, the presence of vertical shear increases the size of the diffusive range while decreasing the slope of the spectrum in most of that range.

We also find that the diffusive range extends to scales substantially larger than the simple estimate *κ*_{e}/*s*, where *s* is a typical stretching rate. This, together with the finite scale of the forcing, implies that there may not be a clear convective range in the stratosphere and might explain why a *k*^{−1} spectrum is not observed. Using observed winds and values of 10^{−1}–10^{−2} m^{2} s^{−1} for the diffusivity corresponding to recent observational estimates (Waugh et al. 1997; Balluch and Haynes 1997; Alisse et al. 2000; Sparling and Bacmeister 2001; Joseph et al. 2003), we derive, even for the lower values of diffusivity, a spectrum that is significantly steeper than *k*^{−1}, even at relatively large scale, and is, in fact, not inconsistent with the observed *k*^{−2} form at scales of 20–200 km.

## 2. Advection–diffusion in the Batchelor regime

### a. Modeling approach

We start with the advection–diffusion equation for the concentration *C*(**x**, *t*) of a passive tracer,

where **v**(**x**, *t*) is the advecting velocity field and *κ* the diffusivity, and we focus on the chaotic-advection regime or Batchelor regime where the velocity field is smooth.

In this regime it may be self-consistently assumed that the typical scales of variation of **v** are much larger than those of *C*; formally,

where ɛ ≪ 1 is a small parameter. It is convenient to work directly with the slow coordinates ɛ**x** and ɛ*t*: we therefore make the substitution (**x**, *t*) → (**x**, *t*)/ɛ in (2.1); we also rescale the diffusivity according to *κ* → ɛ*κ* to obtain

and regard *C* as a function of the slow coordinates.

A useful quantity to characterize the tracer variance is the Wigner function (or action density), defined as

(see Ryzhik et al. 1996; Antonsen et al. 1996). This quantity can be interpreted as the variance of the tracer at a scale given by the wavevector **k**, at location **x** and time *t.* An alternative expression for *ω* is derived using the Fourier transform *Ĉ* of *C,* with

where an asterisk denotes the complex conjugate.

We are now in position to derive an evolution equation for *ω*: differentiating (2.3) with respect to *t* and using (2.2) yields, after some manipulations involving (2.4)–(2.5) and the incompressibility condition ∇·**v** = 0 (Antonsen et al. 1996),

where **v̂** is the Fourier transform of **v**. Using the scale separation assumption ɛ ≪ 1, we find at leading order in ɛ the differential evolution equation

where ∇_{k} is the gradient in **k** space. This equation is central to our analysis; its left-hand side describes the transport of *ω* in the six-dimensional (**x**, **k**) space, while its right-hand side describes the damping of *ω* due to diffusive effects. A solution can be found by integration along characteristics, that is, by integrating the following system of ordinary differential equations:

with initial conditions **x**(0) = **x**_{0}, **k**(0) = **k**_{0}, and *ω*(0) = *ω*_{0}. The first two equations are recognized as the ray equations derived using a WKB approximation in a nondispersive system with phase velocity **v**(**x**, *t*); the third corresponds to the (diffusive) damping of the action.

We solve (2.6) using an ensemble approach: an ensemble of solutions of (2.7) (orbits) is followed in time. Of course, the initial conditions of the ensemble must be chosen to sample the initial distribution of *ω* in (**x**, **k**) space. When required, an approximation to *ω* may be estimated using a binning procedure. This constructs the probability density function (pdf) *p*(*ω, ***x**, **k**, *t*) such that *p*(*ω, ***x**, **k**, *t*) *dωd***x***d***k** is the probability to find tracer with Wigner function in [*ω, **ω* + *dω*] located in the volume [**x**, **x** + *d***x**] × [**k**, **k** + *d***k**] in (**x**, **k**) space, at time *t.* Statistical information about the tracer distribution is then given by the moments of *p*(*ω, ***x**, **k**, *t*). Global measures such as spectra represent averages over **x** and therefore depend only on the **x**-integrated pdf *p*(*ω, ***k**, *t*) = ∫ *p*(*ω, ***x**, **k**, *t*) *d***x**. This is an advantage since it means that a much smaller number of orbits is necessary for accurate numerical calculation than would be the case for quantities depending on details of the full pdf *p*(*ω, ***x**, **k**, *t*).

In practice, we shall always assume that there is no correlation between **x**_{0} and either **k**_{0} or *ω*_{0}. Thus **x**_{0} may simply be regarded as a label for a particular orbit in the ensemble and any space dependence of the pdf will have no significance.

It may be useful to give a physical interpretation of this approach. Each orbit may be interpreted as describing a fluid parcel carrying tracer with a typical wavenumber **k**. As the parcel is advected by the velocity field **v**, it is stretched as a result of the strain field ∇**v**, and the wavevector **k** evolves accordingly. The Wigner function associated with the parcel can be thought of as the local tracer variance and characterizes the amplitude of the tracer fluctuations. The evolution of **x**, **k**, and *ω* is depicted schematically in Fig. 1.

Motivated by application to the stratosphere, we follow Haynes and Anglade (1997) in neglecting the vertical component of the velocity, carrying out the advection in two-dimensions (ideally on an isentropic surface) but taking into account the vertical shear of the horizontal winds. Two spatial coordinates are therefore considered for **x**, while three components are retained for **k**, since its vertical component, *m* say, evolves as a result of the shear. This is seen from (2.7) and the form of the velocity gradient tensor, which is

Although the evolution of the horizontal wavevector **k**_{h} decouples from that of *m,* it is crucial to evaluate *m* since |*m*| ≫ |**k**_{h}| for atmospheric parameters, and thus *m* plays the controlling role in the diffusion process (Haynes and Anglade 1997).

### b. Tracer spectrum

The tracer spectrum is the main focus of this paper. Given the aim of explaining observed horizontal spectra, we consider a tracer spectrum as a function of the modulus of the horizontal wavenumber *k* = |**k**_{h}|, the result of integrating the full spectrum over all orientations of **k**_{h} and over all vertical wavenumbers *m.* There is a relatively straightforward relation between this integrated spectrum and that which would be obtained by averaging over several measurements along horizontal sections (e.g., aircraft flight tracks) in different directions and at different levels. Appendix A discusses this relation in more detail.

To derive the spectrum, we compute the pdf *p*(*ω, **k, **t*), evaluated from the full pdf *p*(*ω, ***x**, **k**_{h}, *m, **t*) by integration over **x**, *m,* and over all orientations of **k**_{h}, and note that the spectrum is the first moment of *p,* regarded as a function of *k* and *t*; that is,

In numerical computations *F̃*(*k, **t*) is calculated directly by summing the values of *ω* for all orbits whose values of *k* fall into a given bin.

The approach described gives the spectrum for an initial value advection–diffusion problem. However, as stated earlier, our principal interest lies in situations where a forcing continuously introduces tracer variance at large scale. The spectrum *F*(*k, **t*) corresponding to the forced situation may be derived from *F̃*(*k, **t*) using the latter as a kind of Green function. Assuming that the forcing is constant in time and that the velocity field is statistically time independent so that *F̃*(*k, **t*) depends only on the time interval *t* since the initial instant, not the initial time itself, yields the spectrum for the forced situation in the form

(Antonsen and Ott 1991). The steady-state spectrum may be obtained by taking the limit *t* → ∞. We will drop the time argument of *F* to denote that this limit has been taken.

The form of the spectrum may be usefully interpreted in terms of an &ldquo=uivalent time” *τ,* defined by

Noting from (2.7) that *ω* = *ω*_{0} exp[−2*κk*(*t*)^{2}*τ*], we interpret *τ* as the time that diffusion would take to reduce the local variance from its initial value *ω*_{0} to its time-*t* value *ω* if the wavenumber had the constant value *k*(*t*). The interest of considering *τ* is that, in stretching flows, it is determined by the recent stretching history along the characteristics, in contrast to *k*(*t*), or equivalently the finite-time stretching rate log[*k*(*t*)]/*t,* which are determined by the entire stretching history. Provided that the correlation time for the stretching is much less than the time required for the exponential decay factor to decrease substantially, it follows (Antonsen et al. 1996) that *τ* becomes independent of *k*(*t*) for large times and, furthermore, that its probability distribution tends to a time-independent form, *M*(*τ*) say. (This holds in most of the cases we investigate later, although not in the Batchelor limit considered in section 3b.) In the two-dimensional (2D) case, it may then be shown that, for large *k,*

where *A* is a constant. The derivation is detailed in appendix B. An expression similar to (2.11) is given by Yuan et al. [2000, their Eq. (16)], but certain details differ and we believe our result to be correct; see appendix B.

This may all be generalized to the layerwise two-dimensional case of interest here by defining, in addition to *τ,* the second equivalent time

The diffusive decay of *ω,* given by *ω* = *ω*_{0} exp[−2*κk*(*t*)^{2}(*τ* + *τ*_{υ})], is controlled by both *τ* and *τ*_{v}. Because the ratio *m*(*t*)/*k*(*t*) is determined by the recent stretching history (Haynes and Anglade 1997), arguments similar to those invoked by Antonsen et al. (1996) suggest that, for large *t,* both *τ* and *τ*_{υ} will be independent of *k.* (Again this does not hold in the Batchelor limit.) In the large-shear limit relevant to the real stratosphere it will be the case that typical values *τ* and *τ*_{υ} of the equivalent times satisfy *τ*_{υ} ≫ *τ.* In the range 1/*κ**τ*_{υ} ≪ *k* ≪ 1/*κ**τ*, the decay factor is well approximated by exp(−2*κk*^{2}*τ*_{;gy}) and the forced tracer spectrum is approximately given by

where *M*_{υ}(*τ*_{υ}) is the pdf of *τ*_{υ}.

In the 2D case, it has been noted that an enhanced variability of the strain leads to a broader distribution of *τ* and hence to a smoother dissipative cutoff (Kraichnan 1968; Yuan et al. 2000). This effect is expected to be much more marked in 3D: while the distribution of *τ* is determined by the Lagrangian statistics of the horizontal strain field, the distribution of *τ*_{υ} is determined by the statistics of both the horizontal strain and the vertical shear. As a result, one might expect a greater spread in the values of *τ*_{υ} than in the values of *τ* and, as a consequence, a dissipative cutoff that is significantly less sharp in the 3D case than in the 2D case. This is confirmed by explicit results obtained below for simple random-strain models. While our discussion will be focussed on the tracer spectrum, we note that it could equally well be formulated in terms of the equivalent-time distributions *M*(*τ*) and *M*_{υ}(*τ*_{υ}) as the Laplace transform formulas (2.11) and (2.12) indicate.

## 3. Random-strain models

Before using realistic stratospheric winds it is interesting to analyze and discuss the behavior of the tracer spectrum using simple models for the velocity field. An important class of such models are random-strain models, which are derived assuming that the strain tensor ∇**v** experienced by a fluid parcel is a random function of time. For the model described in section 2, this implies that the strain tensor is a random function along the characteristics defined by (2.7). Such models have traditionally been applied to turbulent flows (e.g., Kraichnan 1968), but are also relevant to chaotic advection flows with simple spatial structure. Indeed, the chaotic nature of particle trajectories implies that the strain may be quasi-random along Lagrangian trajectories in spite of the Eulerian structure of the flow. Ngan and Shepherd (1999a) have suggested that random-strain models may be particularly relevant to the stratospheric surf zone where the Lagrangian correlation of the strain decreases rapidly with time.

Random-strain models have recently received renewed attention in the turbulence community [e.g., Falkovich et al. (2001) and references therein]. An important result is the form of the spectrum for scales large enough to be unaffected by diffusion (convective range): the stationary spectrum in that range is given by a *k*^{−1} power law (the Batchelor spectrum), irrespective of the specific random process chosen to model the strain (Kraichnan 1968, 1974). This result relies only on the exponential stretching of the wavevectors **k** and on the conservation of *ω*; it can be expected to hold for a wide variety of flows, even if the strain is not random (Pierrehumbert 1992; Antonsen et al. 1996). It also applies (now to the horizontal wavenumber spectrum) when the effect of the vertical shear is taken into account. However, the shape of the spectrum in the diffusive range, that is, the range of scales for which diffusion is important, as well as the extent of this range, depend on the particular strain model and are affected by the vertical shear. These are the aspects that we investigate here.

We follow Haynes and Anglade (1997) and take a strain tensor of the simple form

where *a*(*t*), *b*(*t*), *c*_{1}(*t*), and *c*_{2}(*t*) are independent random processes. This (partly) symmetric form neglects the influence of horizontal rotation as well as possible correlations between the horizontal strain and the vertical shear. We also assume statistical isotropy, taking the same statistics for *a* and *b,* and for *c*_{1} and *c*_{2}, respectively. [Note that the form of (3.1) implies that ∇**v** is not statistically homogeneous; it is alternatively possible to enforce homogeneity, but this imposes a nonzero vorticity fixed by the variance of *a.*]

In the numerical computations presented below, we model *a, **b, **c*_{1}, and *c*_{2} by Ornstein–Uhlenbeck processes (e.g., Gardiner 1985). Each process is then characterized by its variance and its inverse correlation time. We denote by *γ*^{2} the variance of *a* and *b,* and by *σ* the corresponding inverse correlation time; *γ*^{2}_{c} denotes the variance of *c*_{1} and *c*_{2}, and *σ*_{c} the corresponding inverse correlation time. The (Ito) stochastic differential equation obeyed by *a,* for instance, then reads

where *B* is a Brownian process. Correspondingly, the stationary correlation function of *a* is given by

with similar expressions for *b, **c*_{1}, and *c*_{2}.

### a. Kraichnan limit

Analytic results about the tracer spectrum can be obtained when the random processes defining the strain have very short correlation times. In the limit of infinitesimal correlation time, generally called the “Kraichnan limit,” the precise definition of the random processes is unimportant; a single quantity, the time integral of the correlation, suffices to characterize each process. We now derive an equation for the tracer spectrum in this limit. It is obtained by rescaling the strain tensor according to

where 0 < *δ* ≪ 1 and *a, **b, **c*_{1}, and *c*_{2} are independent random functions of their argument satisfying

*η* and *λ* being *O*(1) constants. In (3.2), the correlation times of the various random processes are *O*(*δ*^{2}), and the *O*(*δ*^{−1}) scaling of the strain is such that the corresponding time-integrated correlations are *O*(1). Note that the assumption of a symmetric strain tensor is unimportant here because the rotational part of the strain tensor plays no role in the Kraichnan limit (Chertkov et al. 1995).

In the limit *δ* → 0, Eq. (2.7) for the wavevector **k** simplifies into an equation with white-noise forcing terms. We derive this equation in appendix C and deduce the following elliptic equation for the stationary tracer spectrum *F*(*k, **m*):

This equation is an extension, including the effect of the vertical shear, of that derived by Kraichnan (1968, 1974). In general, it can be solved only numerically; *F*(*k, **m*) may then be integrated with respect to *m* to obtain the horizontal wavenumber spectrum *F*(*k*).

The equation for *F* relevant in the 2D case is obtained by setting *λ* = 0 and *m* = 0 in (3.3). The asymptotic form of its solution for small wavenumbers *k* ≪ *η*/(2*κ*) is *F* ∼ *k*^{−1} (the Batchelor spectrum) and that for large wavenumbers *k* ≫ *η*/(2*κ*) is *F* ∼ *k*^{−1/2} exp(−2*κ**k*/*η*) (Kraichnan 1974).

To analyze (3.3) in the 3D case, with *λ* ≠ 0, it is convenient to use polar coordinates (*r, **θ*) defined by

and to introduce *y*(*r, **θ*) = *kF*(*k, **m*). Equation (3.3) is then rewritten as

where the only remaining parameter is the aspect ratio *α* = *λ*/*η.* This equation can be solved with boundary conditions fixing the value of *y* on a curve *r* = *r*_{0}(*θ*) ≪ 1 (corresponding to forcing at small spatial wavenumbers) and requiring decay for *r* → ∞.

A number of properties of the solution of (3.4) may be deduced analytically. First, for *r* sufficiently small, the solution with the slowest decay in *r* is given by *y* ∼ const; that is, *F* ∼ *k*^{−1}: the Batchelor spectrum is thus recovered for such wavenumbers. Assuming that *α* ≫ 1, as has been argued to be most relevant to the stratosphere, this solution ceases to be a good approximation when *r* ∼ *α*^{−1}, that is, when *k* ∼ *α*^{−1}*η*/(2*κ*). Note that this wavenumber is smaller, by a factor *α*^{−1}, than the diffusive wavenumber in the 2D case. This is a manifestation of the effect of the vertical shear. At this horizontal wavenumber there is a substantial proportion of realizations for which the vertical wavenumber is sufficiently large for vertical diffusion to have a significant dissipative effect (even though the horizontal wavenumber is not large enough for horizontal diffusion to be effective). This remains true as *k* increases until *k* ∼ *η*/(2*κ*) when horizontal diffusion too becomes important.

Asymptotic calculations for the range *α*^{−1}*η*/(2*κ*) ≪ *k* ≪ *η*/(2*κ*), or equivalently *α*^{−1} ≪ *r* ≪ 1, are presented in appendix D and predict the form of the spectrum as

where numerical prefactors have been omitted. This approximation breaks down when *k* ∼ *η*/(2*κ*). Further analytic progress is possible when *k* ≫ *η*/(2*κ*) and details are given in appendix D. It may be shown that

where *ρ*(*α*) = [(1 + 4*α*^{2})^{1/2} − 1]/4. The three asymptotic forms of the spectrum in the Kraichnan limit are summarized in Table 1.

We have verified these asymptotic results by performing numerical experiments using Ornstein–Uhlenbeck processes for the component of the strain tensor. The scaling of the Kraichnan limit is achieved when the variance and correlation times of the processes satisfy *γ*^{2} ≫ 1, *σ* ≫ 1, with *η*^{2} = 2*γ*^{2}/*σ* = *O*(1), and similarly *γ*^{2}_{c} ≫ 1, *σ*_{c} ≫ 1, with *λ*^{2} = 2*γ*^{2}_{c}/*σ*_{c} = *O*(1).

Here, we report results obtained using the values *γ* = 10, *σ* = *σ*_{c} = 200, *γ*_{c} = 200, and a diffusivity *κ* = 1 so that *η* = 1, *λ* = 20, and *α* = 20. This choice of parameters ensures that the range of validity of (3.5), here 0.025 ≪ *k* ≪ 0.5, is wide and well separated from the range of validity of (3.6), here *k* ≫ 0.5. Figure 2 compares (a) the spectrum obtained numerically, (b) the asymptotic result (3.5), and (c) the asymptotic result (3.6). [The prefactors in (3.5)–(3.6) have been adjusted for a good match with the numerical result.] A *k*^{−1} power law is also shown. The figure confirms the validity of the asymptotic estimates; the spectrum is in fact remarkably well approximated for all *k* by a composite of the asymptotic predictions.

The Kraichnan limit provides a good test for the idea that an equivalent diffusivity *κ*_{e} can be employed to mimic the effect of the shear in purely 2D computations. The typical aspect ratio *m*/*k* derived by Haynes and Anglade (1997) is *λ*/*η* = *α,* so the equivalent diffusivity is estimated as

Using this equivalent diffusivity amounts to replacing *m*^{2}/*k*^{2} by *α*^{2} in (3.3). Integration of (3.3) with respect to *m* eliminates the term proportional to ∂^{2}_{mm}(*kF*) and shows that the horizontal spectrum *F*(*k*) is simply the Kraichnan 2D spectrum with *κ*_{e} in place of *κ.* This spectrum is also displayed in Fig. 2: curve (d) shows the spectrum obtained numerically; and curve e shows the large-*k* form, namely,

valid for *k* ≫ *η*/(2*κ*_{e}). The equivalent-diffusivity model predicts correctly the wavenumber at which diffusive effects becomes important, that is, when the spectrum departs significantly from the *k*^{−1} nondiffusive form. However, it is obvious that the spectrum obtained with the equivalent diffusivity does not approximate correctly for all *k* the one derived by taking the shear into account, notably because the diffusive decay is too strong in the former. Indeed the large-*k* form (3.8) clearly disagrees with both (3.5) and (3.6). Note that (3.6) is actually closer to the form of the 2D spectrum with actual, rather than equivalent diffusivity, in the sense that the exponential terms agree, although the algebraically decaying prefactor is much more rapidly decaying in (3.6). It is not difficult to conclude that it is not possible to find a value of *κ*_{e} that mimics, in the 2D spectrum, the effect of the shear in a wide range of *k.*

### b. Batchelor limit

The Batchelor (1959) limit, which assumes very large correlation times for the random processes *a, **b, **c*_{1}, and *c*_{2} in (3.1), is the opposite extreme to the Kraichnan limit. More precisely, the correlation time for the velocity gradient tensor is large compared to both the stretching time scale and the time taken for deformation of the spatial structure from the forcing scale to the dissipation scale. This limit is easily tractable analytically: the slow variation of *a, **b, **c*_{1}, and *c*_{2} allows one to integrate the evolution equation for **k**_{h} and *m* in a deterministic fashion. Quantities such as the power spectrum are thus computed first for each realization of the random processes, then averaged over the statistics of the processes.

The ray equations (2.7) are first rewritten as

where *c* is the projection of (*c*_{1}, *c*_{2}) on **k** (cf. Haynes and Anglade 1997). The large correlation time allows *a, **b,* and *c* to be treated as constants within each realization of the flow. [This shows that in this particular limit the variable *τ* defined by (2.10) and the variables *τ* and *τ*_{υ} cannot be regarded as independent of *k.*]

The deterministic equation for the stationary spectrum *F*(*k, **m*) within each realization, or equivalently for *y* = *kF,* is derived readily. It has the form

and can be solved analytically using characteristics. Defining

we write the solution as

where *h*(·) is an arbitrary function that is determined using the boundary conditions. A simple solution is obtained by taking the boundary condition *y* ∼ *δ*(*m*) for *k* → 0. It follows that

where *C* depends on *a, **b, **c* only.

We now compute the average of *y* over all possible flow realizations, assuming that the variables *a, **b, **c* are normally distributed, that is, taking the probability density function

for *a,* and similarly for *b* and *c.* We further take *γ*_{a} = *γ*_{b} = :*γ.* The choice of *C,* which governs the dependence of the tracer variance forcing on the random processes *a, **b, **c* is important and should be guided by physical considerations. Following Kraichnan (1968) we require that the rate at which the tracer variance is dissipated is independent of *a, **b, **c*; that is, we require that

is independent of *a, **b, **c* and find *C* ∼ (*a*^{2} + *b*^{2})^{−1/4}*c*^{−1} from (3.9). Introducing this result into (3.9) and integrating with respect to *a, **b, **c,* and *m* leads to an expression for the horizontal spectrum of the form

where *α* = *γ*_{c}/*γ* is again the aspect ratio and

As in the Kraichnan limit, there are two asymptotic regimes of interest when *α* ≫ 1, in addition to the *k*^{−1} obtained when diffusion can be neglected. The first corresponds to the range *α*^{−1}*γ*/*κ* ≪ *k* ≪ *γ*/*κ* where diffusion enters only through its effect in the vertical; it leads to the simple power-law spectrum:

[Note that this result depends on the particular (Gaussian) statistics chosen for the random strain and shear.] The other asymptotic regime of interest corresponds to the range *k* ≫ *γ*/*κ* where horizontal diffusion cannot be neglected; the spectrum in this case is given by

Note that this expression makes no assumption on the size of *α.* For *α* = 0 it reduces to the asymptotic form of the 2D spectrum. The three asymptotic forms of *F*(*k*) in the Batchelor limit are summarized in Table 1.

As in the previous section, we confirm our analytical results numerically using Ornstein–Uhlenbeck processes for the strain-tensor components. The validity of the Batchelor limit imposes the constraints *γ, **γ*_{c} ≫ *σ, **σ*_{c} on the variances and inverse correlation times of these processes. In addition, the correlation times are long compared to the time of cascade from the forcing scale *k*_{0} to the diffusive scale *γ*/*κ*, as is required in the Batchelor regime, provided that *σ*^{−1}, *σ*^{−1}_{c} ≫ log(*κ**k*^{2}_{0}/*γ*)/(2*γ*). Parameters satisfying these constraints have been chosen for numerical simulations; they are *γ* = 1, *γ*_{c} = 100, *σ* = *σ*_{c} = 10^{−20}, *k*_{0} = 10^{−5}, and *κ* = 1. Figure 3 shows (a) the spectrum obtained numerically, together with (b) the analytic form (3.10), and (c) the asymptotic result (3.13). The asymptotic power laws are also indicated. The agreement between numerical and theoretical spectra is good. The various asymptotic formulas are seen to provide useful approximations in their range of validity: these are *k* ≪ 0.01 for the *k*^{−1} power law, 0.01 ≪ *k* ≪ 1 for the *k*^{−2} power law, and *k* ≫ 1 for formula (3.13).

Like the Kraichnan limit, the Batchelor limit can be used to assess the concept of equivalent diffusivity. Since the typical aspect ratio between vertical and horizontal scales is *α* as before, (3.7) again provides an estimate for the equivalent diffusivity. The spectrum obtained with this equivalent diffusivity is shown in Fig. 3, curve (d) is derived from numerical calculation, and curve e from the analytical expression for the 2D spectrum, that is (3.10) with *α* = 0 and *κ*_{e} in place of *κ.* As in the Kraichnan limit, the equivalent-diffusivity model predicts approximately the wavenumber at which the spectrum starts to depart significantly from the *k*^{−1} law, but it does not provide a satisfactory approximation for the shape of the spectrum obtained with shear. Again, the large-*k* behavior is markedly different with, according to (3.13), a decay with shear that is essentially the same as the 2D decay but with diffusivity *κ* and not *κ*_{e}.

### c. Finite correlation time

We conclude this section on random-strain models by examining the tracer spectrum when the processes defining the strain and shear have finite correlation times. When this is the case, that is, when *γ*/*σ* and *γ*_{c}/*σ*_{c} are order-one quantities, analytical progress is difficult and we resort to numerical simulations. Theoretical results have been obtained with finite correlation time, which shown in particular that stretching factors (defined by log|**k**|/*t*) follow a large-deviation distribution in the long-time limit (Chertkov et al. 1995). However, the relationship between tracer spectrum and stretching factor is a subtle one, and simple predictions for the spectrum form do not appear to exist.

For the numerical calculations, we have taken parameter values that are roughly representative of the typical values expected in the midlatitude lower stratosphere. This allows comparison between the spectrum obtained here and that obtained in the next section using observed stratospheric winds, thus providing some assessment of the ability of random-strain, Ornstein–Uhlenbeck models to capture the main features of stratospheric mixing. The values chosen are *σ* = *σ*_{c} = 1.15 × 10^{−5} s^{−1}, that is, correlation times of about 3 days as suggested by Ngan and Shepherd (1999b); *γ* = 5.5 × 10^{−6} s^{−1}, that is, an inverse strain of about 2 days; and *γ*_{c} = 1.4 × 10^{−3} s^{−1}. With these parameters, the average (horizontal) stretching rate (or maximum Liapunov exponent)

and the aspect ratio

calculated numerically using a ensemble of 5000 realizations of the Ornstein–Uhlenbeck processes are found to be *h* = 3.5 × 10^{−6} s^{−1} ≈ 0.3 day^{−1} and *α* ≈ 250, respectively. These values are in agreement with those obtained using observed winds in Haynes and Anglade (1997) and in the next section, and they have motivated the choice of *γ* and *γ*_{c}. The aspect ratio *α,* in particular, is roughly given by the ratio of the typical strain to the typical shear *γ*/*γ*_{c}, as expected (Haynes and Anglade 1997). (Note that typical stratospheric values of the strain and shear are likely to be larger than the values of *γ* and *γ*_{c}. This is because vorticity, neglected in our random-strain model but nonzero in the stratosphere, reduces the stretching rate.)

The spectrum is computed using a forcing wavenumber *k*_{0} = 10^{−6} m and a diffusivity *κ* = 10^{−2} m^{2} s^{−1}. Figure 4 displays the spectrum obtained after about 30 days when an approximately stationary state is attained. The spectrum calculated in the absence of shear but with the equivalent diffusivity *κ*_{e} = *κ*(1 + *α*^{2}) = 625 m^{2} s^{−1} is also shown. As in the asymptotic regimes of sections 3a–3b, the spectrum obtained with shear departs from the Batchelor *k*^{−1} spectrum at a wavenumber somewhat smaller than the relevant estimate *h*/*κ*_{e} = 7.5 × 10^{−5} m^{−1}, that is, at a scale somewhat larger than 13 km. This scale is well predicted by the equivalent diffusivity model, although the rest of the spectrum is not. Because of the large value of *α,* a large portion of the spectrum, roughly for *k* for between *h*/*κ*_{e} and *h*/*κ* = 2 × 10^{−2} m^{−1}, is in the regime where diffusion is effective only through its damping of vertical structures.

The conclusion drawn in the asymptotic cases is seen to hold here: the use of an equivalent diffusivity overestimates the diffusive decay of the spectrum, and it is not possible to tune the value of the equivalent diffusivity to reproduce the spectrum obtained with shear satisfactorily. The reason for this behavior, which we have observed for various parameter choices, is simple to understand. The random shear introduces a strong variability in the vertical wavenumber *m,* hence in the diffusion experienced by fluid parcels. This variability means that orbits with a given horizontal wavenumber *k* have their variance distributed over a wider range of values than is the case in the absence of shear. In particular, the number of orbits with large horizontal wave-vector **k** and significant tracer variance is larger than what it is when an equivalent diffusivity is used. The overall result of the shear is therefore a shallower spectrum in the diffusive range, as well as a smoother transition between the *k*^{−1} Batchelor spectrum and the steeper diffusive region. This is also consistent with our asymptotic results (3.6) and (3.13) showing that the effect of the vertical shear on the shape of the large-*k* tail of the spectrum is relatively weak (since this part of the spectrum is dominated by realizations that have experienced little vertical shear).

## 4. Results based on observed winds

The results of the previous section have been obtained using rather arbitrary random processes, chosen for simplicity rather than realism. However, we expect our qualitative conclusions about the tracer spectrum to apply to a wide class of flows. This includes flows for which the strain and shear experienced along fluid-parcel orbits is not random in a strict sense, but simply sufficiently disordered in a loose sense. Stratospheric flows are expected to belong to this category (Ngan and Shepherd 1999a,b). We do not study here the precise applicability of random-strain models to the stratosphere; rather, we concentrate on the effect of stratospheric winds on passive tracers and employ the WKB technique of Antonsen et al. (1996) described in section 2 to compute typical tracer spectra using observed winds. A crude comparison with the results of section 3 allows us to verify that the spectrum obtained is consistent with the idea that the stratospheric strain and shear are random in a Lagrangian sense (this would not be the case if, for instance, Lagrangian trajectories were regular in part of the flow). On the other hand, we can test whether the hypotheses made in section 2 (quasi-two-dimensional flow, scale separation, fixed diffusivity) are relevant to the stratospheric case by comparing the spectrum that we obtained with that derived from aircraft measurements (Strahan and Mahlman 1994; Bacmeister et al. 1996; Sparling and Bacmeister 2001). These measurements suggest that, in the range 20–500 km, the spectrum has a power-law behavior with an exponent close to −2.

As those of Haynes and Anglade (1997), our computations are performed using wind fields derived from stratospheric sounding unit (SSU) geopotential data assuming geostrophic balance. The advection is carried out on the 50-mb isobaric surface rather than on an isentropic surface (with respect to which vertical motion would be weaker), but we are confident that this makes little difference in the spectrum. The strain and shear are computed by finite differences, using winds on the 100-mb and 20-mb surfaces for the latter. Equations (2.7) are integrated for an ensemble of 841 orbits with initial positions **x**_{0} distributed in the latitude band 30°–60°N. The horizontal wavenumber of the forcing is *k*_{0} = 10^{−6} m^{−1} (corresponding to spatial scales of about 6000 km), while the vertical wavenumber is *m*_{0} = 0. Since we are interested in the shape of the spectrum in the diffusive range and in the extent of this range, the choice of the value of the diffusivity *κ* is crucial. This diffusivity must be a turbulent diffusivity, parameterizing the effect of neglected small-scale advection and stretching. For the computations that take into account the vertical shear and whose results are presented in Figs. 5 and 6, we have taken *κ* = 10^{−2} m^{2} s^{−1}. This value is suggested by a number of recent calculations (Waugh et al. 1997; Balluch and Haynes 1997; Alisse et al. 2000; Sparling and Bacmeister 2001; Joseph et al. 2003), and is rather small compared to various other estimates.

Figure 5 shows *kF*(*k*) as a function of *k* (in log–linear coordinates) after 6, 12, 18, 24, and 30 days starting on 10 January 1987. After 30 days, the spectrum has reached a steady state. A striking feature of this steady spectrum is the extent of the diffusive region: there is a significant departure from the *k*^{−1} spectrum (a horizontal line) for wavenumbers as small as 10^{−5} m^{−1} (scales as large as 600 km). Since our results have to be taken with caution for the largest scales for which there is no good separation between velocity and tracer scales, this may explain why a *k*^{−1} power law does not appear in spectra derived from observations. Figure 6 (curve a) displaying *F*(*k*) after 30 days in log–log coordinates confirms that the *k*^{−1} behavior of the spectrum is limited to very large scales. It also shows that, although definitely not a power law, the spectrum is relatively well approximated by a *k*^{−2} law for *k* between 2.5 × 10^{−5} m^{−1} and 5 × 10^{−4} m^{−1}. At smaller scales it is significantly steeper and roughly consistent with that obtained by Sparling and Bacmeister (2001).

It is instructive to compare the spectrum in Fig. 6 with that obtained in section 3c using a random-strain model with realistic parameters and displayed in Fig. 4. Although the two spectra are roughly similar, there are noticeable differences: the departure from the *k*^{−1} power law occurs at larger scale in Fig. 6 than it does in Fig. 4, and the diffusive decay at small scales is shallower. These features cause the spectrum in Fig. 6 to be closer to the observed *k*^{−2} form than the random-strain spectrum of Fig. 4. They also suggest that the processes chosen to represent horizontal shear and vertical strain in the random-strain model do not have a sufficient variability. A reason for this might be the assumption of a symmetric strain tensor (3.1) corresponding to a vanishing vorticity of the flow. For fixed strain and shear, the first effect of a nonzero vorticity is to reduce the horizontal and vertical stretching rates. Our random-strain calculation accounts for this effect in an averaged sense, since we have tuned the values of the strain and shear variance to obtain realistic stretching rate and aspect ratio, but it does not account for the impact of the fluctuations in the vorticity amplitude. This impact is likely to be an extension to large scales of the dissipative range, and a shallowing of the spectrum in this range, consistent with our results.

In section 3, we have emphasized the crucial role played by the vertical shear in shaping the tracer spectrum. Its importance in the stratosphere is also illustrated by curved in Fig. 6, which has been obtained by neglecting the shear and using an equivalent diffusivity. Formula (3.7) has been used with the calculated aspect ratio *α* = 250, yielding *κ*_{e} = 625 m^{2} s^{−1}. Comparison between results obtained with and without shear again shows that the use of an equivalent diffusivity overestimates the dissipative decay at small scales, leading to a spectrum much steeper than the one found when the shear is included explicitly. This somewhat subtle effect of the shear explains why there is no sharp (exponential) dissipative cutoff at small scales in the observed spectra, but rather a moderate steepening (consistent with a *k*^{−2.5} power law; Sparling and Bacmeister 2001). This steepening occurs at a scale comparable to the diffusive scale based on *κ*_{e} and a typical horizontal strain rate of 0.3 days^{−1}, found to be 13 km.

We remark that the results obtained with vertical shear support estimates for the diffusivity *κ* of the order of 10^{−2} m^{2} s^{−1}. Indeed, the spectrum that we obtain with this value is not inconsistent with spectra derived from observations, given the many assumptions made. A value significantly larger for *κ,* by contrast, leads to a spectrum too steep for intermediate scales (*k* ∼ 10^{−5}–10^{−3} m^{−1}) to be consistent with observations; a value significantly smaller leads to a *k*^{−1} law reaching small scales, which is also inconsistent with observations. The latter point is illustrated by Fig. 7 showing the spectrum obtained after 30 days taking into account the vertical shear but with several values for the diffusivity. For the molecular value, *κ* = 10^{−4} m^{2} s^{−1}, the *k*^{−1} spectrum reaches scales ten times smaller than it does for *κ* = 10^{−2} m^{2} s^{−1}. It appears therefore that the Batchelor spectrum is not observed in the stratosphere because the (turbulent) diffusivity is too large for a convective range (where both the influence of the forcing and of the diffusion can be neglected) to exist.

The spectrum *F*(*k*) considered so far is a 2D isotropic spectrum, obtained by integration over all orientations of the horizontal wavevector **k**_{h}. As pointed out earlier, this spectrum differs from the one-dimensional spectrum *G*(*k*) obtained by sampling the tracer field along a horizontal section (as is the case for aircraft measurements). The relationship between *F*(*k*) and *G*(*k*) is discussed in appendix A. Equation (A.2), which allows the computation of *G*(*k*) from *F*(*k*), shows that *G*(*k*) depends on *F*(*l*) only for *l* ≥ *k.* This, and the rapid decay of *F*(*k*) for large *k,* makes it possible to estimate *G*(*k*) in applications where *F*(*k*) is known only for a finite range of *k.* For completeness, we show in Fig. 8 a comparison between the spectrum *F*(*k*) obtained using observed winds (same as in Fig. 6), and the corresponding spectrum *G*(*k*). Here *G*(*k*) is seen to be somewhat steeper than *F*(*k*), particularly at large scale, so that a clear *k*^{−1} power law does not appear even near the forcing scale.

## 5. Discussion

This paper considers the factors that influence the stationary spectrum of stratospheric passive tracers in a forced–dissipative situation. An important feature of the stratospheric flow, that it be dominated by large-scale winds, places the problem in the Batchelor regime of passive advection. This makes it possible to use an ensemble approach based on the WKB approximation (Antonsen et al. 1996) to solve the advection–diffusion equation over a large range of scales, from the planetary scale of the forcing to the diffusive scale. A further simplification is afforded by the quasi-two-dimensional nature of the stratospheric winds, which allows the determination of the horizontal tracer spectrum on a single isentropic surface. Crucially, however, the vertical (or more precisely cross-isentropic) shear of the horizontal wind and its impact on the vertical structure of the tracer are taken into account; this is essential because the vertical contribution dominates the dissipative process (Haynes and Anglade 1997).

The WKB approach provides a computationally efficient technique to calculate the spectrum and, potentially, other statistical diagnostics for a passive tracer that is forced and dissipated. It may be regarded as a complement to the contour-advection and other trajectory-based techniques (which are restricted to conservative, initial value problems but give detailed information about the spatial—as opposed to spectral—distribution of the tracer) and, therefore, as a promising tool to examine how changes in the stratospheric flow, such as seasonal changes, affect tracer distributions.

In the Batchelor regime, the evolution of the tracer statistics depends only on the Lagrangian statistics of the strain tensor. This justifies the use of random-strain models in which the strain components are modeled by simple random processes. With such models, significant analytical progress is possible in the Kraichnan and Batchelor limits corresponding to small and large correlation times of the random processes. These limiting cases are useful to gain some understanding of the influence of the various parameters on the tracer spectrum and, in particular, of the role of the vertical shear. They are not, however, directly relevant to the stratospheric context; we therefore consider a more realistic model, with strain components given by Ornstein–Uhlenbeck processes, that allows the finite correlation times observed in the stratosphere.

The paper goes beyond these idealized models and directly uses analyzed stratospheric winds to compute particle trajectories and the strain components along these trajectories. The comparison between the tracer spectrum so obtained (Fig. 6) and that obtained with a random-strain model with parameters thought to be realistic (Fig. 4) reveals the limitations of our random-strain model. A number of features of that model can be invoked to explain these limitations, among them the use of a partly symmetric strain tensor (assuming an irrotational velocity field) and the neglect of possible correlations between the various components of the strain. More sophisticated random-strain models without these simplifying features could, in principle, be developed. But, because of the difficulties in estimating the many parameters involved, this is arguably an exercise of limited interest. This view is strengthened by the remarkable simplicity with which the WKB approach can be implemented for observed winds.

Several conclusions about stratospheric tracer spectra can be drawn from the different models used in the paper. First, with observed winds, the influence of dissipation extends to scale significantly larger than the diffusive scale *κ*_{e}/*h*, estimated to be of the order of 13 km by Haynes and Anglade (1997) for a diffusivity *κ* = 10^{−2} m^{2} s^{−1}. This can be attributed to the variability in the strain and shear (cf. Yuan et al. 2000). Together with the finite scale of the forcing, this results in a spectrum in which the nondissipative *k*^{−1} behavior is confined to the largest scales. Since these are also the scales for which the WKB approach is invalid, this suggests that the *k*^{−1} power law that appears at large scales in our simulations is not relevant to observed stratospheric spectra.

A second conclusion concerns the effect of the vertical shear on the horizontal tracer spectra. Haynes and Anglade (1997) suggested that this effect may be simply parameterized in purely 2D computations by increasing the diffusivity from its real value *κ* to the equivalent value *κ*_{e} = (1 + *α*^{2}), where *α* is the average aspect ratio of tracer structures. Our results indicate that this parameterization is useful only to estimate the largest scales at which dissipation affects the spectrum so that it departs from the *k*^{−1} form. The diffusive decay predicted using an equivalent diffusivity turns out to be much steeper than the one obtained by taking proper account of the vertical structure of both the tracer and the winds. In fact, there is a large range of scales, a priori between *κ*/*h* and *κ*_{e}/*h* but in practice extending to scales significantly larger than *κ*_{e}/*h*, in which the spectrum, although affected by diffusion, remains relatively shallow. This helps to explain why observed tracer spectra do not exhibit a sharp cutoff at the diffusive scale *κ*_{e}/*h* ∼ 13 km as might have been expected from purely 2D calculations.

Most of our discussions rely on the estimate of the diffusivity *κ* = 10^{−2} m^{2} s^{−1}, which is believed to be realistic. Computations of the tracer spectrum with different values for *κ* lend some support to this idea: when *κ* differs radically from 10^{−2} m^{2} s^{−1}, the spectrum is inconsistent with the observations of a near *k*^{−2} power law. The spectrum computed with *κ* = 10^{−2} m^{2} s^{−1} is roughly comparable to the observed ones in the range of validity of the WKB approach. There is no precise match, however, neither with the spectra of Bacmeister et al. (1996) nor with the arguably more accurate ones of Sparling and Bacmeister (2001). The differences, however, are most striking at small scales, say scales smaller than 20 km, where there is still considerable uncertainty about the representativity of the spectra derived from observations (Sparling and Bacmeister 2001). Furthermore, at such small scales, one might question whether diffusion is a good parameterization of the stratospheric small-scale mixing. As pointed out in Vanneste and Haynes (2000), mixing in the stratosphere is confined within localized, short-lived patches of three-dimensional turbulence and, as a result, highly intermittent and not necessarily well modeled by diffusion. It will be of interest to examine the consequences of this particular type of mixing for stratospheric tracer spectra.

To conclude, we remark on the initial-value problem that describes the evolution of a passive tracer suddenly released in the stratosphere. Antonsen et al. (1996) use the ensemble approach to predict the decay, in a 2D flow, of the tracer variance in such problem. Exactly the same approach could be used for the case studied here, of horizontal flow with vertical shear. However, the usefulness of the prediction for the initial value problem is doubtful. Sukhatme and Pierrehumbert (2002) and Fereday and Haynes (2003, manuscript submitted to *Phys. Fluids*) argue that ensemble approaches for the initial value problem break down at large times and therefore cannot calculate the asymptotic decay. There may be some skill at early times, but only if the initial condition is dominated by spatial scales that are much smaller than the scale of the flow (and homogeneously distributed on the large scale). It is difficult to see this as relevant to the lower stratosphere. What is likely to be the case in the actual initial value problem is that the decay at intermediate times, but not the decay rate at large times, depends on the vertical shear. Application of the ensemble approach to the initial value problem suggests that the intermediate-time decay can be largely captured by the equivalent diffusivity, but whether or not this is actually the case remains to be seen.

## Acknowledgments

The authors thank John Methven, Jai Sukhatme, and Ray Pierrehumbert for comments. This research started while JV was holding an E. C. Marie Curie Fellowship at the University of Cambridge. PHH's research is supported by NERC Grants GR3/12531 and GST/02/2798 and the Isaac Newton Trust; JV's research is supported by a NERC Advanced Fellowship.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

_{2}O and O

_{3}tracer spectra in the stratospheric surf zone.

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

_{2}O measurements, 2, Tracer variability and diabatic meridional circulation.

**,**

**.**

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### 2D Isotropic Spectrum versus 1D Spectrum

In this appendix we compare two definitions of the spectrum as a function of a scalar (horizontal) wavenumber: *F*(*k*) obtained by integration over all orientations of the wavevector **k** and *G*(*k*) obtained in a specific spatial direction (e.g., along an aircraft flight track).

Let *f*(**k**) be the Fourier transform of the concentration correlation; that is,

where the angle brackets denote an ensemble average and spatial homogeneity is assumed. In two dimensions with **k** = *k*(cos*θ,* sin*θ*), *F*(*k*) is defined by

Here *G*(*k*) is the equivalent of *f* for a fixed direction **e**; it is therefore defined by

Comparing *F*(*k*) and *G*(*k*) is meaningful only if the concentration is statistically isotropic. In this case *f*(**k**) = *f*(*k*), *F*(*k*) = *kf*(*k*), and

where *J*_{0}(·) is a Bessel function of order 0. By isotropy, the left-hand side is independent of the direction of **r**; taking **e** for this direction leads to

The 2D isotropic spectrum *F*(*k*) is thus related to the 1D spectrum *G*(*k*) in a nontrivial way: the two functions are integral transforms of each other.

Practical formulas giving *F*(*k*) as a function of *G*(*k*) and vice versa may be derived from (A.1). Inverting the Fourier transform gives

The derivation of the converse expression uses the inversion formula for the Hankel transform (Erdelyi et al. 1954). Its application to (A.1) leads to

After integration by parts the integral with respect to *l,* we can change the order of integration and find

using the Hankel transform of *r*^{−1/2} sin(*lr*) (Erdelyi et al. 1954). Further manipulations yield the simple form

Equations (A.2)–(A.3) emphasize the subtle connection existing between the two power spectra *F*(*k*) and *G*(*k*). Clearly they depend nonlocally on each other: determining *F,* say, at a given scale requires knowledge of *G* at all smaller scales. Note, however, that *F*(*k*) and *G*(*k*) are proportional when they are exact power laws.

### APPENDIX B

#### Derivation of (2.11)

where *R*(*k, **τ, **t*) is the joint pdf for *k* and *τ* at time *t.* For each realization of the flow, **k**(*t*), the wavevector at time *t* is related to its initial value **k**_{0} by a linear map. Taking the maximum eigenvalue of the corresponding matrix to be exp(*ht*) defines the finite-time stretching rate *h.* For *k* ≫ *k*_{0}, we then have

where *ϕ* is the angle between **k**_{0} and the eigenvector associated to the eigenvalue exp(*ht*) (Antonsen et al. 1996). Because the form of the spectrum is insensitive to the details of the distribution of *ϕ,* we assume that *ϕ* is uniformly distributed in (0, 2*π*). We can then rewrite (B.1) in terms of the joint pdf for *h* and *τ* at time *t, **Q*(*h, **τ, **t*) say, as

The corresponding stationary forced spectrum *F*(*k*) is obtained by integration with respect to *t.* Integrating the *δ* function gives

Now, if we assume that, as is generally the case, *Q*(*h, **τ, **t*) decreases rapidly for large *h,* the dominant contribution to the integral comes from moderate values of *h,* typically *h* ∼ *h*, where *h* is the mean stretching rate (or Liapunov exponent). Correspondingly, the large-*t* form of *Q*(*h, **τ, **t*) can be used in (B.2) for wavenumbers *k* ≫ *k*_{0}. As discussed in section 3, this limiting form is particularly simple: for *h**t* ≫ 1, *h* and *τ* become statistically independent, and the pdf of *τ* does not depend on time. Thus, *Q*(*h, **τ, **t*) factorizes according to

The form of *P*(*h, **t*) for *h**t* ≫ 1 is provided by the large-deviation theory (Chertkov et al. 1995; Falkovich et al. 2001). In the range |*h* − *h*| = *O*(*t*^{−1/2}), which gives the dominant contribution to the integral in (B.2), the simple central-limit form, given by

for some constant *a* > 0, can be employed. Introducing (B.3)–(B.4) into (B.2) and integrating with respect to *h* and *ϕ* leads to (2.11), or explicitly

Note that the independence of *h* and *τ* and the central-limit theorem used to derive (B.5) require the additional assumption *σt* ≫ 1, where *σ* is a typical (Lagrangian) inverse correlation time of the strain. This assumption holds for the various models considered in this paper, with the notable exception of the random-strain model in the Batchelor limit examined in section 3b, which has *σ* = 0.

We remark that Yuan et al. (2000) derive an expression similar to but different from (B.5). Their Eq. (16) is formally identical to (B.5) but with *M*(*τ*) replaced by a different function, *M*_{Y}(*τ*) say, which may be written as

where *M*(*τ̂*, *η*) is the joint pdf for the equivalent time *τ* and instantaneous stretching rate *η* = *h*^{−1}*dh*/*dt.* We believe our result (B.5) to be correct. There are difficulties with the derivation of Yuan et al. (2000), possibly related to their assumption of independence between *η* and *k*; consequences of these difficulties are obvious if the strain is such that *M̂*(*τ,* 0) ≠ 0, or if *k* is nondifferentiable, as is the case for random-strain models in the Kraichnan limit considered in section 3a.

### APPENDIX C

#### Spectrum in the Kraichnan Limit

In the short-correlation-time limit *δ* → 0, the components of the strain tensor (3.2) reduce to white-noise processes. The reduced form of (2.7) may be then derived directly from a theorem of Papanicolaou and Kholer (1974) [see also Gardiner (1985), section 6.5]; using the convention of the Ito calculus, we find

where *B*_{1}, *B*_{2}, and *B*_{3} are independent Brownian processes. Replacing *k*_{1 }*dB*_{1} + *k*_{2 }*dB*_{2} by the equivalent process *kdB* (Gardiner 1985, p. 107), leads to a closed equation for *k* &lsqb=uivalent to (A.5) in Haynes and Anglade (1997)]. The complete equations for the tracer Wigner function are then written as

while the Fokker–Planck equation for the associated joint pdf *p*(*ω, **k, **m, **t*) takes the form (Gardiner 1985, p. 96)

### APPENDIX D

#### Asymptotic Results for the Kraichnan Limit

We derive solutions of (3.4) in the limits of large *α* and of large *r.*

First we write *s* = *αr* and take the limit *α* → ∞ in which case (3.4) reduces to

We now seek a solution of the form *y* = exp*G*(*s, **θ*), implying that

This suggests a leading-order solution at large *s* of the form *G*(*s, **θ*) = *sg*(*θ*), where

We seek a solution with *g*′(0) = 0 allowing *g,* and hence *y,* to be an even function of *θ.* There is no analytic solution, but a power series solution for small *θ* may easily be derived as *g*(*θ*) = ±*θ*^{2}/2 + *O*(*θ*^{4}). The negative sign is chosen to give decay away from *θ* = 0. Seeking higher order terms in the expansion for *G*(*s, **θ*) reveals a singular behavior at *θ* = 0. This is a manifestation of the fact that, for small *θ,* neglected terms on the right-hand side of (D.2) become as large as those on the left side and motivates the examination of an inner region where *θ* ∼ *s*^{−1/2}.

Defining a new variable *ξ* = *θs*^{1/2} and substituting *y* = *Y*(*s, **ξ*) into (D.1) implies that

A possible dominant balance involves the three terms on the left-hand side and neglects all the terms on the right-hand side. Solutions may be found by separating variables and, under the assumption of decay as *ξ* → ±∞, may be shown to be of the form

where *n* is a nonnegative integer and He_{n} is the *n*th order Hermite function. Choosing the negative sign in the exponent, to give decay as *s* increases, and focusing on *n* = 0, corresponding to the solution that decays least rapidly (in *s* and in *ξ*), we pose the expansion

with *μ* unknown. Then, at leading order,

implying that *Y*_{0} = *C*_{0} exp(−*ξ*^{2}/2) and at next order,

This equation has bounded solutions for *Y*_{1} only when the right-hand side satisfies the condition

and it follows that *μ* = ½.

We thus have the leading-order solution for *Y* in the inner region as

This, and higher-order terms, may be matched to solutions in the outer region (noting that an expansion in powers of *s*^{−1/2} is needed in that region). However this inner solution contains all the information needed to calculate the leading-order form of the horizontal wavenumber spectrum *F*(*k*) since the major contribution to the integral of *y* with respect to *m* comes from the inner region. Noting that *m* = *s*^{1/2}*ξ* in this region and integrating *Y* with respect to *ξ,* it follows that

Now we consider the regime where *r* ≫ 1. We seek a solution to (3.4) of the form *y* = *z* exp[*rH*(*θ*)], implying that

If we further assume that *z* = *z*_{0}(*θ*)*r*^{μ} + *z*_{1}(*θ*)*r*^{μ−1} + ⋯, then the leading-order term in the equation is *O*(*r*^{μ}) and is given by

As before, we seek a solution that is an even function of *θ,* with *H*′(0) = 0, such that *y* decays as *θ* increases. It is useful to write the solution in the form of a power series in *θ* as *H* = −1 − *ρθ*^{2} + ⋯, where *ρ* = [(1 + 4*α*^{2})^{1/2} − 1]/4. The next term in the equation is *O*(*r*^{μ−1}) and implies that

This equation has a singular point at *θ* = 0 unless *μ* = −*ρ,* in which case, the corresponding equation for *z*_{1}(*θ*) has a singular point at *θ* = 0. The singularity is resolved by an inner region in which *θ* ∼ *r*^{−1/2}. The form of the solution in this inner region is

where *ξ* = *r*^{1/2}*θ, **ν* is an unknown constant to be determined and the leading-order equation is

Applying the boundary condition that *Z*_{0} increases only algebraically as *ξ* → ∞ (which allows matching to the outer region) gives the eigenvalue condition *ν* = −*n*/2 − (2*n* + 1)*ρ,* where *n* is a nonnegative integer and *Z*_{0} ∼ *ξ*^{n} as *η* → ∞. Matching to the outer region gives *μ* = −(2*n* + 1)*ρ.*

Again the largest contribution to the integral of *y* with respect to *m* comes from the neighborhood of *m* = 0 and hence from the neighborhood of *θ* = 0. The integral does not therefore depend on the precise form of the function *z*_{0}(*θ*) and indeed is dominated by the contribution from the inner region. Evaluating the *m* integral by integrating the solution in the inner region with respect to *ξ* and noting that the largest contribution is obtained when the maximum value of *ν,* that is, *ν* = −*ρ,* is chosen, the horizontal wavenumber spectrum follows as

## Footnotes

*Corresponding author address:* J. Vanneste, School of Mathematics, University of Edinburgh, King's Buildings, Mayfield Road, Edinburgh EH9 3JZ, United Kingdom. Email: J.Vanneste@ed.ac.uk