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 m2 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 m2 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 N2O and O3 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 m2 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,

 
tC + v · ∇C = κ2C,
(2.1)

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,

 
v = vx, ɛt) and C = C(x, t),

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

 
tC + v · ∇C = ɛ2κ2C
(2.2)

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

 
formula

(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

 
formula

A dual expression to (2.3) is obtained by substituting (2.4) into (2.3) to find

 
formula

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),

 
formula

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

 
formula

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:

 
formula

with initial conditions x(0) = x0, k(0) = k0, 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ωdxdk is the probability to find tracer with Wigner function in [ω, ω + ] located in the volume [x, x + dx] × [k, k + dk] 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) dx. 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 x0 and either k0 or ω0. Thus x0 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.

Fig. 1.

Interpretation of the orbit equation (2.7): a parcel of fluid starting at x0 with tracer at scale k0 and Wigner function ω0 is advected at time t to a new position x where the tracer scale is k and the Wigner function ω. The reduction in amplitude of ω is indicated by the change from black to gray

Fig. 1.

Interpretation of the orbit equation (2.7): a parcel of fluid starting at x0 with tracer at scale k0 and Wigner function ω0 is advected at time t to a new position x where the tracer scale is k and the Wigner function ω. The reduction in amplitude of ω is indicated by the change from black to gray

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

 
formula

Although the evolution of the horizontal wavevector kh decouples from that of m, it is crucial to evaluate m since |m| ≫ |kh| 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 = |kh|, the result of integrating the full spectrum over all orientations of kh 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, kh, m, t) by integration over x, m, and over all orientations of kh, and note that the spectrum is the first moment of p, regarded as a function of k and t; that is,

 
formula

In numerical computations (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 (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 (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

 
formula

(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

 
formula

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,

 
formula

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

 
formula

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κk2τ;gy) and the forced tracer spectrum is approximately given by

 
formula

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

 
formula

where a(t), b(t), c1(t), and c2(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 c1 and c2, 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, c1, and c2 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; γ2c denotes the variance of c1 and c2, and σc the corresponding inverse correlation time. The (Ito) stochastic differential equation obeyed by a, for instance, then reads

 
da = −σa dt + γ(2σ)1/2 dB,

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

 
a(t)a(0)〉 = γ2 exp(−σ|t|),

with similar expressions for b, c1, and c2.

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

 
formula

where 0 < δ ≪ 1 and a, b, c1, and c2 are independent random functions of their argument satisfying

 
formula

η 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):

 
formula

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 Fk−1 (the Batchelor spectrum) and that for large wavenumbers kη/(2κ) is Fk−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

 
formula

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

 
formula

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 = r0(θ) ≪ 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, Fk−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 α−1r ≪ 1, are presented in appendix D and predict the form of the spectrum as

 
F(k) ∼ exp[−2(2ακk/η)1/2],
(3.5)

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

 
F(k) ∼ k−1/2−ρ(α) exp(−2κk/η),
(3.6)

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

Table 1.

Asymptotic forms of the tracer spectrum F (k) in the Kraichnan and Batchelor limits for α ≫ 1. See text for the definitions of the various parameters. Note that the results in the Batchelor limit apply only to Gaussian strain and shear statistics

Asymptotic forms of the tracer spectrum F (k) in the Kraichnan and Batchelor limits for α ≫ 1. See text for the definitions of the various parameters. Note that the results in the Batchelor limit apply only to Gaussian strain and shear statistics
Asymptotic forms of the tracer spectrum F (k) in the Kraichnan and Batchelor limits for α ≫ 1. See text for the definitions of the various parameters. Note that the results in the Batchelor limit apply only to Gaussian strain and shear statistics

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 γ2c ≫ 1, σc ≫ 1, with λ2 = 2γ2c/σ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.

Fig. 2.

Log–log plot of the horizontal spectrum in the Kraichnan limit with a vertical shear: (a) spectrum obtained numerically (see text for the parameter values); (b) asymptotic form (3.5); (c) asymptotic form (3.6). The spectrum obtained by modeling the vertical shear using an equivalent diffusivity is also shown: (d) spectrum obtained numerically; (e) asymptotic form (3.8). A k−1 power law is indicated

Fig. 2.

Log–log plot of the horizontal spectrum in the Kraichnan limit with a vertical shear: (a) spectrum obtained numerically (see text for the parameter values); (b) asymptotic form (3.5); (c) asymptotic form (3.6). The spectrum obtained by modeling the vertical shear using an equivalent diffusivity is also shown: (d) spectrum obtained numerically; (e) asymptotic form (3.8). A k−1 power law is indicated

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

 
κe = κ(1 + α2).
(3.7)

Using this equivalent diffusivity amounts to replacing m2/k2 by α2 in (3.3). Integration of (3.3) with respect to m eliminates the term proportional to 2mm(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,

 
formula

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, c1, and c2 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, c1, and c2 allows one to integrate the evolution equation for kh 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

 
formula

where c is the projection of (c1, c2) 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

 
formula

and can be solved analytically using characteristics. Defining

 
formula

we write the solution as

 
formula

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

 
formula

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

 
formula

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

 
formula

is independent of a, b, c and find C ∼ (a2 + b2)−1/4c−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

 
formula

where α = γc/γ is again the aspect ratio and

 
formula

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:

 
F(k) ∼ k−2.
(3.12)

[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

 
formula

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 k0 to the diffusive scale γ/κ, as is required in the Batchelor regime, provided that σ−1, σ−1c ≫ log(κk20/γ)/(2γ). Parameters satisfying these constraints have been chosen for numerical simulations; they are γ = 1, γc = 100, σ = σc = 10−20, k0 = 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).

Fig. 3.

Log–log plot of the horizontal spectrum in Batchelor limit: (a) spectrum obtained numerically (see text for the parameter values); (b) analytic result (3.11); (c) asymptotic form (3.13). The spectrum obtained by modeling the vertical shear using an equivalent diffusivity is also shown: (d) numerical result; (e) analytic result. The k−1 and k−2 power laws are indicated

Fig. 3.

Log–log plot of the horizontal spectrum in Batchelor limit: (a) spectrum obtained numerically (see text for the parameter values); (b) analytic result (3.11); (c) asymptotic form (3.13). The spectrum obtained by modeling the vertical shear using an equivalent diffusivity is also shown: (d) numerical result; (e) analytic result. The k−1 and k−2 power laws are indicated

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)

 
formula

and the aspect ratio

 
α = 〈|m|/|k|〉

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 k0 = 10−6 m and a diffusivity κ = 10−2 m2 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 m2 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.

Fig. 4.

Log–log plot of the horizontal spectrum obtained using Ornstein–Uhlenbeck processes with stratospheric parameter values: (a) with vertical shear; (b) modeling the vertical shear by an equivalent diffusivity. The k−1 and k−2 power laws are also shown. The caret indicates the scale of the variance source

Fig. 4.

Log–log plot of the horizontal spectrum obtained using Ornstein–Uhlenbeck processes with stratospheric parameter values: (a) with vertical shear; (b) modeling the vertical shear by an equivalent diffusivity. The k−1 and k−2 power laws are also shown. The caret indicates the scale of the variance source

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 x0 distributed in the latitude band 30°–60°N. The horizontal wavenumber of the forcing is k0 = 10−6 m−1 (corresponding to spatial scales of about 6000 km), while the vertical wavenumber is m0 = 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 m2 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.

Fig. 5.

Log–linear plot of kF(k) vs k after 6, 12, 18, 24, and 30 days obtained using stratospheric data for the wind strain and shear, and a diffusivity κ = 10−2 m2 s−1. The caret indicates the scale of the variance source

Fig. 5.

Log–linear plot of kF(k) vs k after 6, 12, 18, 24, and 30 days obtained using stratospheric data for the wind strain and shear, and a diffusivity κ = 10−2 m2 s−1. The caret indicates the scale of the variance source

Fig. 6.

Log–log plot of F(k) vs k after 30 days: (a) with shear and for a diffusivity κ = 10−2 m2 s−1; (b) without shear and for an equivalent diffusivity κe = 625 m2 s−1. The k−1 and k−2 power laws are also shown

Fig. 6.

Log–log plot of F(k) vs k after 30 days: (a) with shear and for a diffusivity κ = 10−2 m2 s−1; (b) without shear and for an equivalent diffusivity κe = 625 m2 s−1. The k−1 and k−2 power laws are also shown

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 m2 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 m2 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 m2 s−1, the k−1 spectrum reaches scales ten times smaller than it does for κ = 10−2 m2 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.

Fig. 7.

As in Fig. 6, but with diffusivity κ = (a) 10−4 m2 s−1, (b) 10−3 m2 s−1; (c) 10−2 m2 s−1, and (d) 10−1 m2 s−1

Fig. 7.

As in Fig. 6, but with diffusivity κ = (a) 10−4 m2 s−1, (b) 10−3 m2 s−1; (c) 10−2 m2 s−1, and (d) 10−1 m2 s−1

The spectrum F(k) considered so far is a 2D isotropic spectrum, obtained by integration over all orientations of the horizontal wavevector kh. 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 lk. 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.

Fig. 8.

Comparison between (a) 1D spectrum G(k) and (b) 2D isotropic spectrum F(k) obtained after 30 days with κ = 10−2 m2 s−1

Fig. 8.

Comparison between (a) 1D spectrum G(k) and (b) 2D isotropic spectrum F(k) obtained after 30 days with κ = 10−2 m2 s−1

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 m2 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 m2 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 m2 s−1, the spectrum is inconsistent with the observations of a near k−2 power law. The spectrum computed with κ = 10−2 m2 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

REFERENCES
Alisse
,
J-R.
,
P. H.
Haynes
,
J.
Vanneste
, and
C.
Sidi
,
2000
:
Estimation of lower-stratospheric diffusivity from microstructure measurements.
Geophys. Res. Lett.
,
27
,
2621
2624
.
Antonsen
,
T. M.
, and
E.
Ott
,
1991
:
Multifractal power spectra of passive tracers convected by chaotic fluid flows.
Phys. Rev.
,
44A
,
851
857
.
Antonsen
,
T. M.
,
Z.
Fan
,
E.
Ott
, and
E.
Garcia-Lopez
,
1996
:
The role of chaotic orbits in the determination of power spectra of passive scalar.
Phys. Fluids
,
8
,
3094
3104
.
Bacmeister
,
J. T.
,
S. D.
Eckermann
,
P. A.
Newman
,
L.
Lait
,
K. R.
Chan
,
M.
Loewenstein
,
M. H.
Proffitt
, and
B. L.
Gary
,
1996
:
Stratospheric horizontal wavenumber spectra of winds, potential temperature, and atmospheric tracers observed by high-altitude aircraft.
J. Geophys. Res.
,
101
,
9441
9470
.
Balluch
,
M.
, and
P. H.
Haynes
,
1997
:
Quantification of lower stratospheric mixing processes using aircraft data.
J. Geophys. Res.
,
102
,
23487
23504
.
Bartello
,
P.
,
2000
:
Using low-resolution winds to deduce fine structure in tracers.
Atmos.–Ocean
,
38
,
303
320
.
Batchelor
,
G. K.
,
1959
:
Small-scale variation of convected quantities like temperature in turbulent fluid. Part 1. General discussion and the case of small convectivity.
J. Fluid Mech.
,
5
,
113
133
.
Chen
,
P.
,
1994
:
The permeability of the Antarctic vortex edge.
J. Geophys. Res.
,
99
,
20563
20571
.
Chertkov
,
M.
,
G.
Falkovich
,
I.
Kolokolov
, and
I.
Lebedev
,
1995
:
Statistics of a passive scalar advected by a large-scale two-dimensional velocity field: Analytic solution.
Phys. Rev.
,
51E
,
5609
5627
.
Erdelyi
,
A.
,
W.
Magnus
,
F.
Oberhettinger
, and
F. G.
Tricomi
,
1954
:
Tables of Integral Transforms.
McGraw-Hill, 391 pp
.
Falkovich
,
G.
,
K.
Gawȩdzki
, and
M.
Vergassola
,
2001
:
Particles and fields in fluid turbulence.
Rev. Mod. Phys.
,
73
,
913
975
.
Gardiner
,
C. W.
,
1985
:
Handbook of Stochastic Methods. 2d ed.
Springer-Verlag, 442 pp
.
Haynes
,
P. H.
, and
J.
Anglade
,
1997
:
The vertical-scale cascade in atmospheric tracers due to large-scale differential advection.
J. Atmos. Sci.
,
54
,
1121
1136
.
Joseph
,
B.
,
B.
Legras
, and
F.
Lefèvre
,
2003
:
Vertical diffusivity in the lower stratosphere from Lagrangian back-trajectory reconstructions of ozone profiles.
J. Geophys. Res., 108, 4562, doi:10.1029/2002JD003045
.
Kraichnan
,
R. H.
,
1968
:
Small-scale structure of a scalar field convected by turbulence.
Phys. Fluids
,
11
,
945
953
.
Kraichnan
,
R. H.
,
1974
:
Convection of a passive scalar by a quasi-uniform random straining field.
J. Fluid Mech.
,
64
,
737
762
.
Methven
,
J.
, and
B.
Hoskins
,
1999
:
The advection of high-resolution tracers by low-resolution winds.
J. Atmos. Sci.
,
56
,
3262
3285
.
Ngan
,
K.
, and
T. G.
Shepherd
,
1997
:
Comments on some recent measurements of anomalously steep N2O and O3 tracer spectra in the stratospheric surf zone.
J. Geophys. Res.
,
102
,
24001
24004
.
Ngan
,
K.
, and
T. G.
Shepherd
,
1999a
:
A closer look at chaotic advection in the stratosphere. Part I: Geometric structure.
J. Atmos. Sci.
,
56
,
4134
4152
.
Ngan
,
K.
, and
T. G.
Shepherd
,
1999b
:
A closer look at chaotic advection in the stratosphere. Part II: Statistical diagnostics.
J. Atmos. Sci.
,
56
,
4153
4166
.
Papanicolaou
,
G.
, and
W.
Kohler
,
1974
:
Asymptotic theory of mixing stochastic differential equations.
Commun. Pure Appl. Math.
,
27
,
641
668
.
Pierrehumbert
,
R. T.
,
1991
:
Large-scale horizontal mixing in planetary atmospheres.
Phys. Fluids
,
3A
,
1250
1260
.
Pierrehumbert
,
R. T.
,
1992
:
Spectra of tracer distribution: A geometric approach.
Nonlinear Phenomena in Atmospheric and Oceanic Sciences, G. F. Carnevale and R. T. Pierrehumbert, Eds., Springer-Verlag, 27–46
.
Pierrehumbert
,
R. T.
,
1994
:
Tracer microstructure in the large-eddy dominated regime.
Chaos, Solitons Fractals
,
4
,
1091
1110
.
Ryzhik
,
L. R.
,
G.
Papanicolaou
, and
J. B.
Keller
,
1996
:
Transport equations for elastic and other waves in random media.
Wave Motion
,
24
,
327
370
.
Shepherd
,
T.
,
J. N.
Koshyk
, and
K.
Ngan
,
2000
:
On the nature of large-scale mixing in the stratosphere and mesosphere.
J. Geophys. Res.
,
105D
,
12433
12446
.
Shraiman
,
B. I.
, and
E. D.
Siggia
,
2000
:
Scalar turbulence.
Nature
,
405
,
639
646
.
Sparling
,
L. C.
, and
J. T.
Bacmeister
,
2001
:
Scale dependence of tracer microstructure: Pdfs, intermittency and the dissipation scale.
Geophys. Res. Lett.
,
28
,
2823
2826
.
Strahan
,
S. E.
, and
J. D.
Mahlman
,
1994
:
Evaluation of the SKYHI general circulation model using aircraft N2O measurements, 2, Tracer variability and diabatic meridional circulation.
J. Geophys. Res.
,
99D
,
10319
10332
.
Sukhatme
,
J.
, and
R. T.
Pierrehumbert
,
2002
:
Decay of passive scalars under the action of single scale smooth velocity fields in bounded two-dimensional domains: From non-self-similar probability distribution functions to self-similar eigenmodes.
Phys. Rev.
,
66E
.
Paper 056302
.
Vanneste
,
J.
, and
P. H.
Haynes
,
2000
:
Intermittent mixing in strongly stratified geophysical fluids as a random walk.
J. Fluid Mech.
,
411
,
165
185
.
Waugh
,
D. W.
, and
R. A.
Plumb
,
1994
:
Contour advection with surgery: A technique for investigating finescale structure in tracer transport.
J. Atmos. Sci.
,
51
,
530
540
.
Waugh
,
D. W.
, and
Coauthors
,
1994
:
Transport of material out of the stratospheric Arctic vortex by Rossby wave breaking.
J. Geophys. Res.
,
99
,
1071
1088
.
Waugh
,
D. W.
, and
Coauthors
,
1997
:
Mixing of polar air into middle latitudes as revealed by tracer-tracer scatter plots.
J. Geophys. Res.
,
102
,
13119
13134
.
Yuan
,
G. C.
,
K.
Nam
,
T. M.
Antonsen
,
E.
Ott
, and
P. N.
Gudzar
,
2000
:
Power spectrum of passive scalars in two dimensional chaotic flows.
Chaos
,
10
,
39
49
.

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,

 
formula

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

 
formula

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

 
formula

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

 
formula

where J0(·) 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

 
formula

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

 
formula

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

 
formula

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

 
formula

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

 
formula

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)

It follows from (2.8) and from (2.10) that the transient spectrum may be expressed as

 
formula

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 k0 by a linear map. Taking the maximum eigenvalue of the corresponding matrix to be exp(ht) defines the finite-time stretching rate h. For kk0, we then have

 
kk0|cosϕ| exp(ht),

where ϕ is the angle between k0 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

 
formula

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

 
formula

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 hh, 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 kk0. As discussed in section 3, this limiting form is particularly simple: for ht ≫ 1, h and τ become statistically independent, and the pdf of τ does not depend on time. Thus, Q(h, τ, t) factorizes according to

 
Q(h, τ, t) = P(h, t)M(τ).
(B.3)

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

 
formula

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

 
formula

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, MY(τ) say, which may be written as

 
formula

where M(τ̂, η) is the joint pdf for the equivalent time τ and instantaneous stretching rate η = h−1dh/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 (τ, 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

 
formula

where B1, B2, and B3 are independent Brownian processes. Replacing k1 dB1 + k2 dB2 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

 
formula

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

 
formula

Multiplying this equation by ω and integrating with respect to ω (including an integration by parts), then with respect to t as required by (2.8), leads to Eq. (3.3) for the stationary spectrum F(k).

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

 
formula

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

 
formula

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

 
g2 + g2 − tan2θ = O(s−1).
(D.2)

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 ξ = θs1/2 and substituting y = Y(s, ξ) into (D.1) implies that

 
formula

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

 
formula

where n is a nonnegative integer and Hen is the nth 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

 
formula

with μ unknown. Then, at leading order,

 
Y0 + (1 − ξ 2)Y0 = 0,

implying that Y0 = C0 exp(−ξ2/2) and at next order,

 
Y1 + (1 − ξ 2)Y1 = ξY0 + (2μ − 1/2)Y0.

This equation has bounded solutions for Y1 only when the right-hand side satisfies the condition

 
formula

and it follows that μ = ½.

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

 
Ys1/2 exp(2s1/2) exp(−ξ 2/2).

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 = s1/2ξ in this region and integrating Y with respect to ξ, it follows that

 
F(k) ∼ exp[2(2ακk/η)1/2].

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

 
formula

If we further assume that z = z0(θ)rμ + z1(θ)rμ−1 + ⋯, then the leading-order term in the equation is O(rμ) and is given by

 
H2 = H2 = 1 + α2 tan2θ.

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

 
formula

This equation has a singular point at θ = 0 unless μ = −ρ, in which case, the corresponding equation for z1(θ) 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

 
z = Z0(ξ )rν + Z1(ξ )rν−1 + ⋯,

where ξ = r1/2θ, ν is an unknown constant to be determined and the leading-order equation is

 
Z0 − (4ρ + 1)ξ Z0 − 2(ρ + ν)Z0 = 0.

Applying the boundary condition that Z0 increases only algebraically as ξ → ∞ (which allows matching to the outer region) gives the eigenvalue condition ν = −n/2 − (2n + 1)ρ, where n is a nonnegative integer and Z0ξn as η → ∞. Matching to the outer region gives μ = −(2n + 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 z0(θ) 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

 
F(k) ∼ k−1/2−ρ(α) exp(−2κk/η).

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