Turbulence-induced fluctuations in ionization and application to PMSE



R. J. Hill1, D. E. Gibson-Wilde2, J. A. Werne2, D. C. Fritts2

1Environmental Technology Laboratory, Environmental Research Laboratories, National Atmospheric and Oceanic Administration, 325 Broadway, Boulder CO 80303, U.S.A.

2Colorado Research Associates, 3380 Mitchell Lane, Boulder CO 80301, U.S.A.





The temporal evolution of a turbulent layer is calculated in detail by solving the hydrodynamic equations. The turbulence is initiated by a Kelvin-Helmholtz instability. The field of potential-temperature fluctuations serves as a tracer for modeling entrainment of the mixing ratios of ionized constituents hypothesized to be present in the upper polar mesosphere. This entrainment modeling provides the input to a turbulence advection model capable of calculating the spectra and cospectra of ions and electrons. The turbulence advection model is used as a sub-grid-scale model and is required because, given present or foreseeable computer capabilities, numerical solutions cannot span the enormous range of spatial scales from the depth of the shear layer to the smallest scales on which the most massive ions diffuse. The power spectrum of electron number-density fluctuations obtained from the turbulence advection model is compared with that measured by a rocket during the STATE (Structure and Atmospheric Turbulence Environment) experiment; agreement is found for a case of massive ions. The radar cross section for Bragg scattering is calculated from the electron number-density power spectrum and is used to calculate the signal-to-noise ratio (S/N) for the Poker Flat 50 MHz radar. The resultant S/N is then compared with the radar measurements obtained during the STATE experiment. These comparisons support the hypothesis that massive ions can cause polar mesosphere summer echoes from turbulent layers. Large-scale morphology of the turbulent layer obtained from rocket and radar measurements is reproduced by the hydrodynamic solution.



1. Introduction

It is desirable to solve the hydrodynamic equations for a given mesospheric dynamical situation to obtain the fields of velocity and of the mixing ratios of entrained mesospheric constituents. Resultant fields can be used to predict radar observables and realizations of measurements along rocket trajectories. As impressive as present-day computational facilities may be, they are insufficient to include both the large-scale atmospheric dynamics that initiate and sustain turbulence and the smallest scales where dissipation smooths the most slowly diffusing constituents. To simulate a Kelvin-Helmholtz billow, we dimensionalize the hydrodynamics solution to the size of the large-scale atmospheric dynamics. Consequently, the Bragg wavelength of VHF radars is smaller than or nearly equal to the grid size used in the hydrodynamics solution and the smallest scale of the most slowly diffusing constituent is much smaller than the grid scale. Therefore, a sub-grid-scale model is needed. Our basic plan is to use numerical solutions to calculate realizations of the hydrodynamic fields using practical grid resolution, and to use those fields as input to a turbulence advection model that provides sub-grid-scale statistics. From the inertial range of the numerical solution to beyond dissipation-range scales, the turbulence advection model calculates the power spectra and cospectra of mesospheric ion mixing ratios. For simplicity we refer to those spectra and cospectra as ion spectra. The ion spectra are statistics (i.e., expectation values), not realizations. The power spectrum of electron number-density fluctuations (referred to as the electron spectrum for simplicity) is obtained from the ion spectra by use of charge neutrality. Evaluation of this electron spectrum at the Bragg wave number of a given radar allows evaluation of the radar scattering cross section. The electron spectrum is compared with that obtained from a rocket measurement of electron number density.

The relevant physical attributes of the charged species are their charge, mass, and momentum-transfer collision frequency; the latter depends on charge, mass, and size. Traditional nomenclature, such as charged aerosol, charged dust, and cluster ion, that differs for different degree of hydration or chemical constituents distracts from the simplicity of the physics that we present. For simplicity we refer to all charged species as ions, in the sense that an ion is a group of atoms having a charge, no matter how large that group may be.

The turbulence advection model employed here is the same as that developed by Hill and Bowhill (1976) and refined by Hill (1978b). There have been many subsequent experimental and theoretical studies that have verified the accuracy of the model and extended its capabilities; these studies are reviewed by Hill and Mitton (1998). The turbulence advection model was first used by Hill and Bowhill (1976) to estimate equatorial radar backscatter power on the basis of large hydrated ions (masses less than 200 nucleons) that cause Schmidt numbers as large as 1.4. Kelley and Ulwick (1988) hypothesized that very massive ions having very large Schmidt numbers (i.e., much larger than 1.4) cause the observed high-wave-number enhancement of the electron spectrum and the radar backscatter power observed during PMSE.

Our present focus is on the phenomenon of PMSE caused by turbulence; some PMSE are observed to be not directly caused by turbulence (Lubken et al., 1993, 1998). Much has been written about PMSE; we refer the reader to recent review articles (Cho and Rottger, 1997; Rottger, 1993, 1994; Cho and Kelley, 1993). Our present hydrodynamics solution is that of a Kelvin-Helmholtz instability producing a billow and subsequently, a turbulent layer. In future studies, we intend to consider breaking gravity waves. The STATE measurements of PMSE provide convenient data on rocket electron number-density profiles and power spectra, data on the dynamics (such as shear, potential-temperature gradients, energy-dissipation rates), and radar signal-to-noise ratio (S/N) from the Poker Flat 50 MHz radar. In Sec. 8, we compare our hydrodynamics solutions and electron spectra with the STATE data. A brief introduction to STATE data is given in Sec. 8.

2. Simulation of Turbulence Caused by Kelvin-Helmholtz Instability

The simulation of Kelvin-Helmholtz (KH) instability was performed using a spectral code that solves the incompressible Navier-Stokes equations in a stream function/vorticity representation. Horizontal boundary conditions in the streamwise direction (along the mean wind direction) and spanwise direction (horizontal and perpendicular to the mean wind direction) are assumed to be periodic. A cosine expansion in the vertical direction and free-slip boundaries at the top and bottom of the computational domain are employed to represent the horizontal velocity components; this enables specification of a velocity profile having zero vertical velocity at the boundaries.

The initial velocity was taken to have vanishing vertical and spanwise components and the streamwise component, V(z), was taken to have the profile

V(z) = U0 tanh(z/h) , (1)

where z is height, h is the shear depth, and U0 is the magnitude of the streamwise component of velocity far above and below the shear layer. The maximum initial shear exists at z = 0 and is U0/h. The initial buoyancy frequency, N, is independent of height. The Richardson, Reynolds, and Prandtl numbers are defined, respectively, by

Ri = N 2/(U0/h)2 , (2)

Re = U0h/ , (3)

and Pr = /D , (4)

where is kinematic viscosity and D is the thermal diffusivity. If Ri, Re, and Pr are assigned values, then the numerical solution can proceed from the two initial conditions, which are N equals a constant and Eq. (1), and from an initial small forcing of vorticity structure. The values assigned for the solution used in our present study are Ri = 0.05, Re = 2000, and Pr = 1. The resulting hydrodynamic solution applies to any values of U0, h, and N that satisfy Eq. (2); for the present solution this is N 2/(U0/h)2 = 0.05. Thus, we are free to assign values to U0 and h to simulate the KH instability in many different flow cases. Choosing values for U0 and h determines N from Eq. (2), and from Eq. (3); then Eq. (4) determines D from .

The KH simulation was performed in a computational domain of dimensions L0, L0/3, and 2L0 in the streamwise, spanwise, and vertical directions, respectively, where L0 is the streamwise domain size. The spanwise extent is adequate to allow vigorous and locally isotopic vortex dynamics at small scales of motion. Sensitivity studies were performed to ensure that the upper and lower boundaries are sufficiently far from the shear layer so as to not influence the shear layer dynamics. The shear depth, h, and spanwise dimension L0 are assigned the ratio L0/h = 4; this is the most numerically efficient way to activate the most rapidly growing eigenfunction according to linear viscous theory (K. Julien, personal communication, 1997). The simulation was initiated at t = 0 with the most unstable two-dimensional (2D) linear-eigenfunction vorticity structure with a maximum amplitude equal to 0.07U0/h. A three-dimensional (3D) noise spectrum in velocity having a maximum vertical vorticity equal to 0.014U0/h was used to seed the transition to 3D structure; the results are insensitive to the details of this noise spectrum. Model resolution varied throughout the simulation as required to resolve: 1) the quasi-2D mean and KH structures, 2) the inertial range character of turbulence at intermediate scales, and 3) several decades of spectral amplitude within the viscous range. During the most vigorous turbulence, the maximum resolution reaches 720, 240, and 1440 grid points in the streamwise, spanwise, and vertical directions, respectively.

In Sec. 8 we use the values U0 = 18.5 m s-1 and h = 270 m. We mention these values here to give quantitative meaning to our variables. For example, Ri = 0.05 then requires that N = 0.0153s-1; the buoyancy period is Tb = 2/N = 410 s, and the billow turn-over time is T0 L0/U0 = 4Ri1/2/N = 2Ri1/2Tb = 0.45Tb = 183 s. Also, L0 = 4h = 3393 m such that the computational domain is 3393 m by 1136 m by 6786 m in the streamwise, spanwise, and vertical directions, respectively. At maximum resolution the grid spacing is L0/720 = 4.71 m in all three directions. For our case study in Sec. 8, at the height of 88.75 km, the kinematic viscosity of air, air, is 1.87 m2 s-1. Given Re =2000, U0 = 18.5 m s-1, and h = 270 m, Eq. (3) gives = 2.5 m2 s-1. It is too expensive to change Re and run the simulation for any given choice of U0 and h such that equals air for a given temperature and pressure. It is also unnecessary because the numerical solution is intended to simulate the large-scale dynamics. Thus, that = 2.5 m2 s-1 whereas air = 1.87 m2 s-1 is of no concern to us, and would remain of no concern even if the difference between and air was much greater than it is. Also of no concern are the facts that the Prandtl number of air is 0.72 whereas the Prandtl number in Eq. (4) has been assigned the value unity and that, for our values of U0 and h, Eq. (4) gives D = = 2.5 m2 s-1, whereas for our case study in Sec. 8 the thermal diffusivity of air at 88.75 km has the value Dair = air/0.72 = 2.6 m2 s-1. Dimensionless time, t, is marked on our first two figures; t is time (in seconds) multiplied by U0/h. Therefore, t = 100 corresponds to 1459 seconds for our choice of U0 = 18.5 m s-1 and h = 270 m.

The fact that for realistic values of U0 and h we have nearly equal to air, (or that we could cause = air by use of either a slightly different Re or U0 and h) is remarkable. It means that our computer simulations replicate the KH dynamics from the billow scale to the viscous scale for altitudes near the mesopause given present day computer resources. In nondimensional terms this means that Re =2000 is a realistic Reynolds number for the KH instability at altitudes near the mesopause.

Cross sections in the streamwise-vertical plane (centered in the spanwise direction) of the vorticity magnitude and potential-temperature perturbation (about a constant mean) are shown at various times spanning the first half of our KH simulation in Figs. 1 and 2. The evolution in Figs. 1 and 2 is seen to be approximately 2D during times less than t 49 (see the top left panel in Figs. 1 and 2 which corresponds to times 2Tb 4T0) At these times there are increasingly sharp potential-temperature and velocity gradients caused by fluid entrainment and baroclinic torques acting on mass-density gradients within the entraining fluid. Thereafter, streamwise-aligned, counter-rotating vortices quickly arise within the outer billow (Klaassen and Peltier, 1985), where the fluid is both convectively and inertially unstable (Fritts et al., 1996; Palmer et al., 1996). The close proximity of the convectively unstable layers and adjacent sheared stable layers leads to orthogonal vortex alignments, loop vortices, and the vortex interactions and twist-wave dynamics that have been seen in related wave-breaking dynamics to drive the turbulence cascade towards ever smaller scales of motion (Arendt et al., 1997; Andreassen et al., 1998; Fritts et al., 1998). The result shown in Fig. 1 is small-scale turbulence that: 1) exhibits a high degree of spatial and temporal intermittency, 2) is confined early on ( t 66 to 83) to the outer billow, 3) achieves the smallest scales and the largest vorticity (and hence largest energy-dissipation rates as well) in the billow core somewhat later ( t = 103 to 146), and 4) thereafter (t > 146) spreads into a more horizontally homogeneous layer having increasing dissipation scales and decreasing energies. The evolution in the potential temperature field in Fig. 2 reveals: 1) the initial 2D banding accompanying billow roll-up (see t = 49), 2) generation of sharp potential-temperature gradients near the top and bottom of the turbulent layer where strong entrainment occurs, 3) small secondary billows along these interfaces constitute continued instability, and 4) strong and intermittent entrainment events at small scales at the interfaces. As the layer restratifies ( t 200), the potential temperature gradients lessen, but remain greatest at the top and bottom of the layer. Werne and Fritts (1999) present other results from the numerical simulation, including spectra of potential temperature and the 3 velocity components, and profiles and temporal evolution of statistics.

3. Multipolar Diffusion

We use the term multipolar diffusion to refer to multiconstituent plasma diffusion. Multipolar is a generalization of the nomenclature ambipolar, the latter is traditionally used to refer to plasma diffusion for the case of electrons with one species of ion.

Multipolar diffusion was studied by Hill and Bowhill (1976,1977) and Hill (1978a). Here, we present a simplified derivation of multipolar diffusion theory. In particular, we do not present the complications caused by nonzero Debye length, nor the nonlinearity of the equations, nor the coupling to potential-temperature fluctuations caused, for instance, by the temperature dependence of ion diffusion coefficients (those complications are considered by Hill and Mitton, 1998). Production and loss of ionization are neglected as are the effects of the geomagnetic and gravitational fields. To simplify the presentation, our notation list is as follows:

is the number density of species .

N is the number density of neutral gas.

is the mixing ratio of species .

P is atmospheric pressure.

T is potential temperature.

is velocity of the neutral gas.

is the velocity of species relative to the neutral gas, wherein is the velocity of species .

is the flux of species relative to the neutral gas.

is the total electric field.

is the charge of species .

e is the elementary unit of charge, i.e., the magnitude of the electron's charge.

is the charge number (including the sign) of species .

is the mass of species .

is the momentum-transfer collision rate for species .

is the mobility of species .



is the diffusion coefficient of the number density of species , wherein TB is absolute temperature multiplied by Boltzmann's constant.

is the current density relative to the flow of neutral gas.

Our equations are (Hill and Mitton, 1998):

(5)



(6)



(7)

Equation (5) expresses the balance of collisional drag on the left-hand side and the electric force and diffusion on the right-hand side. Equation (6) is the continuity equation for the mixing ratios. It is important to obtain equations for the mixing ratios because the numerical simulation can treat mixing ratios as conserved quantities. The electric field caused by multipolar diffusion is given in terms of the net multipolar charge density by Eq. (7).

To specialize the equations to diffusion phenomena, we multiply Eq. (5) by and sum over , which results in an equation for the total current relative to the flow of neutral gas:



(8a)



In plasma diffusion, the electric field retards the diffusion of some species and enhances the motion of other species such that little charge develops and little current flows. Neglecting relative to the charge flow of individual species (i.e., relative to at least some terms in the sum in Eq. (8a)) yields an equation for the electric field produced by the multipolar diffusion, namely (8b)

This equation is valid even if the plasma is devoid of electrons. However, if electrons are sufficiently plentiful to affect charge neutrality, then the electron's term dominates in both sums in Eq. (8b) because the electron mobility, e, is orders of magnitude greater than ion mobilities. That is, Eq. (8b) becomes ; solving for the electric field gives



(8c)



This is essentially Hill's (1978a) equation (22). One can recognize Eq. (8c) as being Eq. (5) applied to the electrons with the electron's collisional drag term (left-hand side of Eq. (5)) neglected. That is, the electrons diffuse so as to enforce approximate charge neutrality. This fact is verified by the yet more general treatment of plasma modes by Hill and Bowhill (1977).

Electron bite-outs have been observed during some PMSE (Pedersen et al., 1970, Ulwick et al. 1988, Inhester et al., 1990, Klostermeyer, 1996, Lubken, 1998). Perhaps the approximation leading to Eq. (8c) does not apply to bite-outs.

Now, Qe is obtained from approximate charge neutrality, that is, setting the left-hand side of Eq. (7) to zero and dividing by the number density of the neutral gas, i.e.,

(9)

The more complete treatment of substituting Eq. (8c) into Eq. (7) results in the Debye shielding formula which differs from Eq. (9) at scales less than or on the order of the Debye length (Hill, 1978a). We substitute Eq. (9) into Eq. (8c) and Eq. (8c) into Eq. (5) to obtain a formula for ion flux, , in terms of the ion mixing ratios, Q . The divergence of this is substituted into Eq. (6) to obtain the continuity equation for the ion mixing ratios. However, the electric force term in Eq. (5) is nonlinear because Eq. (7) contains number densities. Linearizing the resulting continuity equation and denoting fluctuations by primes, the ion continuity equation is

(10)



This continuity equation is used in the next section to derive equations for ion spectra. Because of the linearization, Eq. (10) and the following theory applies only for Q much smaller than <Q> for each .

4. Turbulence Advection Model for Ionization

Here, we use the continuity equation for the ion mixing ratio, Eq. (10), to derive the coupled equations for spectra and cospectra of ion mixing ratios and apply the turbulence advection model. Let {Q} denote the spatial Fourier transform of Q. The Fourier transform of Eq. (10) gives the equation for {Q}. Multiplying the resultant equation by {Q}*, where the superscript * indicates the complex conjugate, and denotes another ion species, one can obtain the equation for the quantity

= RE [{Q}* {Q}],

where RE indicates the real part. The statistical cospectrum of ion mixing ratios (k) is the average of over a locally-isotropic ensemble integrated over spherical shells in wave-vector space, i.e.,

(k) = 4k2 ,

where the angle brackets denote the aforementioned average, and k is the spatial wave number. Averaging the equation for one obtains (Hill and Mitton, 1998) the equation for (k):







(11a)

where

(11b)

and (k) is the spectral transfer function that arises from the term in Eq. (10) that contains the velocity. The relationship between mixing-ratio covariance and (k) is

(12) where dk is the differential solid angle in k-space. Thus, (k) is normalized such that

(13)

Of course, if denotes the same ion species as , i.e., = , then (k) is the power spectrum and the left-hand side of Eq. (13) is the variance of the ion mixing ratio. The term in Eq. (11a) that contains (D + D) is the free-diffusion term, and the term that contains c and c is the coupling term caused by the multipolar electric field.

The presence of two unknown functions, (k) and (k), in Eq. (11a) constitutes the turbulence closure problem and necessitates the use of the turbulence advection model to express

(k) in terms of (k). The turbulence advection model by Hill and Bowhill (1976) is used; namely

(14)



where is energy-dissipation rate per unit mass of fluid. The parameters of this model are determined by comparison with experiments by Hill (1978b), and have been subsequently further confirmed (Hill and Mitton, 1998); the parameters are: = 0.72 (the Obukhov-Corrsin constant), a = 1.4, and k* = 0.074, wherein,

= (air3/)1/4

is the Kolmogorov microscale. For locally stationary turbulence, we neglect the time-derivative term in Eq. (11a). Substituting Eq. (14) into Eq. (11a) gives

(15)





If there are n species of ions, then Eq. (15) yields n(n + 1)/2 coupled, homogeneous, linear differential equations. They are readily solved by using a predictor-corrector algorithm. The boundary conditions in k for solution of Eq. (15) are the inertial-convective-range formulas

(16)



where is the rate of dissipation of the covariance in Eq. (13). Of course, if denotes the same ion species as , i.e., = , then is the rate of dissipation of the variance of mixing-ratio fluctuations of ion species . Values of dissipation rates and are to be obtained from the KH simulation. A value of k within the inertial-convective range is selected; the (k) are assigned the values according to Eq. (16), and the coupled Eqs. (15) are integrated toward higher wave numbers.

In Sec. 7 we need two other types of spectra. One is the power spectrum of potential-temperature fluctuations denoted by TT(k), and the other is the cospectrum of ion mixing ratio and potential temperature, T(k). The equations to calculate TT(k) and T(k) using the turbulence advection model are analogous to Eq. (15); they are obtained from the ion mixing-ratio continuity equation, Eq. (10), and the continuity equation for potential temperature by the same methods used to obtain Eqs. (11a) and (15). Details of that derivation are given by Hill and Mitton (1998). The boundary conditions in k for calculating TT(k) and T(k) are the same as Eq. (16), with , replaced by TT and T, respectively, where TT is the dissipation rate of potential-temperature variance, and T is that of the covariance of ion mixing ratio and potential temperature. The predictor-corrector algorithm must now solve n(n +1)/2 + n +1 equations if there are n species of ions.

5. Schmidt Numbers

Before Eq. (15) for (k) is solved numerically, it is rendered dimensionless by division by and by use of and air. As in Sec. 2, we make the distinction between the kinematic viscosity of air, air, at the temperature and air density of interest and the parameter in Eq. (3). Similarly, the corresponding equations for TT(k) and T(k) are rendered dimensionless. Examining Eqs. (15) and (11b) we see that the ratio D /air must appear in the dimensionless equations. Thus, the Schmidt numbers Sc = air/D naturally appear in our equations. There is one Schmidt number in our coupled equations for each species of ion. There is no Schmidt number for the electrons in our formulation, nor is there an effective electron diffusivity. The electrons distribute themselves as required by approximate charge neutrality because they are by far the most mobile species. An effective electron diffusivity can apply only if all ion species have the same diffusion coefficient; that restrictive case does not apply here. (The most general case for which an effective electron diffusivity applies, including the effect of nonzero Debye length and negative ions, is given in equation (9.5) by Hill and Bowhill, 1976.)

The ion diffusion coefficients are calculated from the momentum-transfer collision frequency of ions colliding with neutral air molecules. The collisions are modeled as the polarization interaction with a rigid-sphere repulsion. The polarization interaction is attractive and of long range, whereas the rigid-sphere repulsion acts at short ranges when the ion and neutral molecule are in contact. This collision model was derived in detail by Langevin (1905) (a modernized translation from the French text is in an appendix by McDaniel, 1964). Langevin's model is implemented for an ion consisting of an ice sphere surrounding a central core and having an arbitrary charge; the ion collides with O2, N2, and Ar. Following Klostermeyer (1996), the mass density of ice is 0.93 g cm-3. Although Langevin's (1905) momentum-transfer collision model is not perfect for mesospheric ions, it is the best model available for use without undue labor.

If the ion's radius is small enough or its charge is large enough, then the polarization interaction dominates over the rigid-sphere repulsion, and vice versa. Ion diffusion coefficients vary from being dominated by the polarization interaction to being dominated by the rigid-sphere repulsion as the number of attached water molecules increases, i.e., as the radius of the ion increases. For the case of a very large ion, the formula for the kinematic viscosity of air divided by the asymptotic formula for the rigid-sphere repulsion gives the following simple formula for the Schmidt number:

Sc = 6.0 nm-2 [2/( + 1)] r2 ,

where r is the radius of the ion , and is absolute temperature divided by 110.4 K. Whereas this asymptotic formula is useful for estimation purposes, we use Langevin's (1905) model to calculate Sc.

6. Coupling the Turbulence Advection Model to the Hydrodynamic Solution

The dissipation rates , , T, and TT are determined from the numerical solution of the hydrodynamic equations; a wave number k in the inertial-convective range is selected; at that point, all values in Eq. (16) are known such that Eq. (16) determines the boundary values for the solution toward higher wave numbers of Eq. (15) and the corresponding equations for TT(k) and T(k).

The dissipation rates are the result of large-scale production of energy and mixing-ratio variance. However, if the Reynolds number is large enough for an inertial range to exist, then the dissipation rates can be calculated by means of the balance of production and dissipation, that is, by employing the direct-dissipation formulas. The numerical simulation described in Sec. 2 does have an inertial range (Werne and Fritts, 1999). The energy-dissipation rate is therefore calculated from the velocity field of the numerical simulation by means of the direct-dissipation formula:

(17)



where the over bar denotes an average over some chosen volume and repeated indices are summed. The dissipation rate of the variance of potential-temperature, TT, is obtained from the direct-dissipation formula:

(18)



where T is the fluctuation of potential temperature calculated from the numerical simulation. Formulas analogous to Eq. (18) are used for calculating and T (Hill and Mitton, 1998). In Eqs. (17) and (18) the values of kinematic viscosity, , and thermal diffusivity, D, are the values used in the numerical simulation as described in Sec. 2. As mentioned in Sec. 2, the values air and Dair that correspond to the case study in Sec. 8 differ somewhat from the values of and D. Nevertheless, Eqs. (17) and (18) yield correct values of dissipation rates for the geophysical dynamics being simulated because altering the values of and D results in altered values of the simulation's gradients such that dissipation rates obtained from Eq. (17) and (18) are unchanged. This property follows from the balance of production and dissipation, which, in turn, follows from the fact that the dissipation rates are caused by large-scale flow properties and that the Reynolds number of the simulation is large enough so that an inertial range exists.

7. Electron Number-Density Spectrum and Radar Cross Section

One objective is to calculate the electron number-density power spectrum (electron spectrum), (k), which is denoted by no subscripts to distinguish it from the power spectrum of electron mixing-ratio fluctuations, ee(k). At this point we have the cospectra and power spectra of the ion mixing ratios. By Fourier transformation of charge-neutrality, Eq. (9), we can immediately obtain

(19a)



Using Eq. (9) we can also obtain the cospectrum of the electron mixing ratio and potential temperature, eT(k),

(19b)

From the definition of mixing ratio, electron mixing-ratio fluctuations can be caused by electron number-density fluctuations and by fluctuations of neutral air number density. Because the Mach number of the KH-initiated turbulence is very small, pressure fluctuations have a negligible effect on air density (Lumley and Panofsky, 1964). Therefore, the ideal gas law requires that the air density fluctuations are inversely proportional to fluctuations of either temperature or potential temperature. To first order in fluctuations we have

(20)



where N, Ne, and T are ambient values of number density of air, electron number density, and potential temperature, and the prime denotes the fluctuations of the quantity. Obtaining the spectrum from Eq. (20) yields

(21)

Now, Eq. (21) gives the desired electron spectrum by use of Eqs. (19a,b) and TT(k).

We digress to illustrate Eq. (21) in Fig. 3 using one of the cases from Sec. 8 (the other cases would appear qualitatively alike if shown in Fig. 3). For the case considered, Fig. 3 shows the 3 terms in Eq. (21), their sum, (k), as well as the Bragg wave number of the Poker Flat radar, Kolmogorov's dissipation wave number

kK = -1 = (/air3)1/4 ,

and the 3 Batchelor's dissipation wave numbers

kB = (/air D2)1/4 = -1Sc1/2

corresponding to the 3 ion Schmidt numbers. Batchelor's dissipation wave number for the temperature fluctuations (i.e., kB = -1Pr1/2, where Pr = 0.72) is not shown because it is nearly the same as kK. All of the spectra have the k-5/3 power law at the left side of Fig. 3. The dotted curve in Fig. 3 has a power law slightly steeper than k-1 in what would otherwise be a viscous-convective range because of the coupling to the multipolar electric field in the last term in Eq. (15). Slightly to the left of the right-most triangle one sees a bulge toward high wave number in the solid curve; this bulge is caused by the heaviest ion in the case considered. In Fig. 3 we see that neither the third nor second terms on the right-hand side of Eq. (21) are significant. Even if they were significant at low wave numbers (which requires much reduced ion dissipation rates relative to those in Sec. 8), they would be negligible at the Bragg wave number indicated on Fig. 3 because eT(k) and TT(k) decrease rapidly in their dissipation ranges relative to ee(k) if at least one ion Schmidt number is very large. Thus, if large ions are present, then only the first term in Eq. (21) contributes significantly to (k) for k kK . In the literature on PMSE, one often finds the level of (k) caused by very large Schmidt numbers compared with what (k) would be if neutral-air fluctuations produced the electron number-density fluctuations, presumably with Qe held constant. The right-most term in Eq. (21) is the level of (k) caused exclusively by neutral-air fluctuations for constant Qe.

To compare our result with rocket data we have another calculation to perform. Within a limited height range, rocket data is essentially data along a straight line. The power spectrum of such data is called the 1D spectrum (k1), where k1 is the wave-vector component along the line. For locally isotropic turbulence, the relationship of (k1) to the 3D spectrum (k) is (Tennekes and Lumley, 1972)



(k1) = k1 k-1 (k) dk. (22)

Another objective is to calculate the radar cross section. Let denote the backscatter cross section per unit volume and per unit solid angle for Bragg scattering from isotropic fluctuations. This cross section is given by:



= re2 22 K-2 (K) , (23)



where re is the classical electron radius, and K is the Bragg wave number, which is related to the radar's frequency f in MHz by K = (4/c) f ×106, wherein c is the speed of light. The cross section in Eq. (23) is the same as was defined by Villars and Weisskopf (1955) and Ottersten (1969). The cross section defined by Royrvick and Smith (1984) is Eq. (23) multiplied by 4. Note that reflectivity is 4.

8. Comparison with STATE Data

The STATE experiment included measurements by the 50 MHz Poker Flat radar, and rocket measurements of electron number density, temperature, and winds (Ulwick et al., 1988; Fritts et al., 1988). The first rocket of the STATE campaign was launched on June 15, 1983 at 0230 UT and was referred to by Ulwick et al. (1988) as STATE 1; it yielded high-quality electron number-density data, but no winds or temperatures. The second (STATE 2) rocket salvo was 3 hrs. 20 min. after the STATE 1 rocket; STATE 2 yielded winds and temperatures but no useable electron number-density data. The STATE 2 data was used by Fritts et al. (1988) to describe the dynamical state of the mesosphere and lower thermosphere. The earlier dynamical situation during the STATE 1 rocket salvo might have been different. The mesopause temperature minimum of 130 K was located near 88 to 90 km, where a maximum in wind shear exceeding 55 m s-1 km-1 was observed (Fritts et al., 1988). Particles as large as those in noctilucent clouds have been observed below the mesopause during such cold mesopause temperatures (Lubken et al., 1996). As the first application of our model, we consider the turbulent layer observed in the altitude range 88 to 89.5 km by the STATE 1 rocket. Kelley and Ulwick (1988) discuss the evidence that this layer is turbulent; we accept their judgement that it was a turbulent layer. Certainly some PMSE are not directly caused by turbulence (Lubken et al., 1993, 1998). In the altitude range 88 to 89.5 km, the STATE 1 rocket trajectory was about 9 km north-north-east of the center of the radar beam (Fig. 1, Ulwick et al., 1988). This turbulent layer produced very large radar backscattered power (Ulwick et al., 1988; Kelley and Ulwick, 1988).

Using STATE data, Watkins et al. (1988) obtained the energy-dissipation rate, , from the radar's spectral width using the radar's 15-degree off-zenith beam direction. The spectral width method is fraught with uncertainties (Hocking, 1983, 1996). None of their values pertain to the specific turbulent layer that we study here. Nevertheless, we mention their values here for qualitative comparison with the simulation's values. Their 4- to 12-hour averaged values were typically less than 0.1 W kg-1. However, they also observed bursts of in layers a few kilometers thick having duration no longer than a few minutes; these bursts included values of as large as 1 W kg-1. A value of in such a burst is more relevant to our investigation of the observed turbulent layer than are averages over many hours.

We selected a single time, t = 146, in the KH simulation to compare with the STATE data. The vorticity and potential-temperature fluctuations for t = 146 are shown in Figs. 1 and 2. Based on typical KH observations, we dimensionalized the KH simulation using U0 = 18.5 m s-1 and h = 270 m. We will retain those values despite the fact that slightly different values would give a better fit to the STATE data because this helps illustrate the uncertainties that we face when comparing with data. One can now compare the greatest initial shear (at the middle height) of the simulation, i.e., U0/h = 68.5 m s-1 km-1, with the shear of about 60 m s-1 km-1 observed by the STATE 2 rocket salvo. We can now determine the buoyancy frequency from Eq. (2) and Ri = 0.05; we obtain N = 0.0153 s-1. Using the relationship of N to the potential-temperature gradient, GT, and using the temperature of 130 K and the gravitational acceleration appropriate to 88 km, we obtain the initial value GT = 3.2 K km-1. The initial gradient of absolute temperature is then -6.3 K km-1. In the height range of the turbulent layer, Fritts et al. (1988) show absolute-temperature gradients from 40 K km-1 to -20 K km-1, although the relevance of the greatly varying temperature gradients measured during the STATE 2 rocket to the gradients that might have existed 3 hrs 20 min earlier during the STATE 1 salvo is unknown.

To examine the vertical morphology of the simulation's potential-temperature field, we obtain a vertical profile of potential temperature from the center of the simulation. To make it relevant to STATE we use a temperature of 130 K at the simulation's middle height, equate temperature and potential temperature at the middle height (i.e., the reference pressure in the definition of potential temperature is taken to be equal to the pressure at the middle of the turbulent layer), use the aforementioned gradient of 3.2 K km-1 to scale the simulation's potential-temperature gradient, and place the simulation's dimensioned turbulent layer at the height of the layer observed by the STATE 1 rocket. The result is the simulation's dimensioned potential-temperature profile shown in Fig. 4. We compare with the electron number-density profile observed by the STATE 1 rocket by tracing that profile from figure 12 of Ulwick et al. (1988) onto our Fig. 4. The potential-temperature profile and the electron number-density profile have similar characteristics. Both have steep gradients at the top and bottom of the turbulent layer. Muschinski and Wode (1998) observe similar steep gradients in temperature and specific humidity at the top and bottom of a turbulent layer in the troposphere.

The similarity of profiles in Fig. 4 does not prove that the STATE 1 turbulent layer was initiated by a KH instability, because a breaking gravity wave would cause a similar morphology after the turbulence had become fully developed. The early stages of instability are when the signature of KH dynamics is most distinct from breaking gravity-wave dynamics. In Fig. 4 the depth of the layer is about 15% greater in potential temperature than in electron number density, which suggests that we should reduce h by about 15% to better fit (at t = 146) the electron number-density profile. At the top and bottom of Fig. 4, the potential-temperature profile retains the initial potential-temperature gradient of 3.2 K km-1. The potential-temperature profile in Fig. 4 is relatively flat at the middle heights of 88.1-89.2 km; in comparison, the electron number-density profile has more variability. Reference to the panel for t = 146 in Fig. 2 shows that, depending on the horizontal position at which one obtains the potential-temperature profile, one could encounter more variability near the top and bottom of the layer than is seen in the potential temperature in Fig. 4.

We have averaged the simulation's energy-dissipation rate and the potential-temperature dissipation rate in rectangular parallelepipeds having depth of 100 m and horizontal extent of 1.13 by 1.13 km. This horizontal extent approximately corresponds to the Poker Flat radar's full two-way beam width multiplied by the range to the layer center of 88.75 km, whereas the 100-m depth is about 1/3 of the radar's range resolution. In Figs. 5 and 6 we show profiles of these averaged dissipation rates. The profile called "center" has the parallelepipeds centered on the vertical line at the streamwise and spanwise middle of the computation domain, whereas the profile called "edge" has this vertical line at the streamwise edge and spanwise middle of the domain (because the horizontal boundary conditions are periodic, a volume at a streamwise edge would be depicted in Figs. 1 and 2 as having half of the volume on the right-hand side of a panel and half on the left-hand side of a panel). The potential-temperature dissipation rates in Fig. 5 have two pronounced maxima at the top and bottom of the turbulent layer. These maxima are caused by entrainment at the strong potential-temperature gradient at the top and bottom of the layer shown in Fig. 4. At the center of the computational domain at t = 146, the potential-temperature dissipation rate is greater at the top of the turbulent layer than at the bottom, whereas at the edge of the domain it is greater at the bottom of the layer. The energy-dissipation rate profiles in Fig. 6 have their maxima in mid layer and much smaller values at the top and bottom of the layer. The values are greater at the center of the computation domain than at the at the edge. The maximum values of 0.12 and 0.4 W kg-1 are within the range of the values observed by Watkins et al. (1988).

For our comparison with STATE data, we use the potential temperature as a tracer for the ion mixing ratios. Since the potential-temperature profile is initially linear with height, we let the ion mixing-ratios be initially linear with height across the vertical extent of the turbulent layer. Any conserved scalar quantities that have proportional initial profiles will have proportional dissipation rates. Thus, the ion mixing-ratio dissipation rates are proportional to that of potential temperature; the coefficient of proportionality is GG /(GT)2, where G is the initial mixing-ratio gradient of ion . That is, = TT GG /(GT)2, similarly, T = TT G /GT. For this simple case we need not solve the mixing-ratio continuity equations during the simulation's computation. On the basis of the data by Ulwick et al. (1988), at 88.75 km the electron number density and its scale height were taken to be 8×103 cm-3 and 4.5 km, respectively; these values are used to require that the electron and ion profiles collectively satisfy charge neutrality.

We tried ion models having both "light" and "heavy" positive ion species. By trial and error we found ion models that produce agreement with the STATE electron power spectrum obtained for the height interval 88 - 89 km. Two ion models are given in Table 1, wherein the scale heights H and the N /Ne satisfy charge neutrality. The gradients are expressed in terms of scales heights which are defined by H = -Q(z0) [G(z0)]-1, where z0 = 88.75 km is the middle height of the layer (this definition does not imply an exponential profile). The ions' radii in Table 1 are based on the mass density of ice being 0.93 gr cm-3 (Klostermeyer, 1996). Ions larger than those in Table 1, apparently of both positive and negative charge, have been observed in the summer polar mesosphere (Havnes et. al., 1996a,b). A motivation for hydration number 15 for the light positive ion case is the measurement of positive-ion hydration number up to 20 by Bjorn et. al. (1985). However, an even lighter positive ion would not change the results from the light positive-ion case.

We now have all that is needed to use the turbulence advection model to calculate (k1) from Eq. (22). In Fig. 7 we show 4 calculations of (k1) obtained using the center and edge of the computational domain and using the heavy- and light-ion models. The dissipation rates were averaged over the height range 87.75-89.33 km as well as over the 1.13 by 1.13 km horizontal area. The model spectra in Fig. 7 approach the k-5/3 inertial-convective range power law at the lowest wave numbers and are, because of the multipolar electric field, slightly steeper than the k-1 viscous-convective range power law at intermediate wave numbers. In Fig. 8 we show the STATE 1 electron number-density 1-D spectrum transcribed from figure 6a by Ulwick et al. (1988); their spectrum is from the height interval 88-89 km. Figure 8 also contains the advection model's (k1) from the light-positive-ion model and the edge of the domain; that is, the long-dashed curves in Figs. 7 and 8 are the same (k1). To avoid too many curves on Fig. 8, only that one curve from Fig. 7 appears in Fig. 8. Our model (k1) lies about a factor of 2 below the STATE 1 spectrum at all wave numbers in Fig. 8. In making this comparison we note that arithmetic averages of the STATE 1 spectral level shown on the logarithmic ordinate in Fig. 8 would appear near the top of the STATE 1 spectral data. The factor of 2 is a consequence of our choice of U0 = 18.5 m s-1 and h = 270 m and the time t = 146 during the simulation. Had we chosen h = 230 m, which is about 15% less than 270 m as previously mentioned for a better fit in Fig. 4, and chosen U0 = 19.8 m s-1, then our model (k1) would be a factor of 2 greater in Fig 8. However, one could also use an earlier time in the simulation when the dissipation rates are greater than the values in Figs. 5 and 6. The rocket trajectory is a single line through the turbulence that might have pierced a region having high or low dissipation rates or some intermediate values. Averaging dissipation rates along vertical lines through the simulation could produce many (k1), some of which might have the same spectral level as the STATE 1 spectral data. The above discussion illustrates uncertainties that we have when comparing with the data.

In their figure 12, Ulwick et al. (1988) show the electron number-density 1-D spectral level evaluated at the Bragg wave number of the Poker Flat radar, namely, (K) in our notation (see Eqs. (22) and (23)). They obtained their (K) by spectral analysis of the STATE 1 rocket-measured electron number-density in many height intervals. In their figure 12, the (K) values are consistent with their subsequent calculation of S/N and with their spectrum in our Fig. 8 only if the (unstated) units of those (K) are (m/rad)(cm-3)2. We change those units to cm-5 and copy their (K) from their figure 12 onto our Figs. 9a,b where we also show profiles of our model's (K), as calculated from the dissipation rates in Figs. 5 and 6. Their (K) have maxima closer to the center of the turbulent layer than do our (K). Perhaps this is a consequence of the particular line along which the rocket pierced the turbulent layer. Our (K) have maxima at almost the same location as does the potential-temperature dissipation rate in Fig. 5, which is also where our have their maxima. In Figs. 9a,b we also show (K), which is proportional to radar cross section, Eq. (23). By differentiating Eq. (22) we have (K) = -n (K), where the so-called spectral index is -n = d(log (k1))/d(log k1) evaluated at k1 = K. The distances between corresponding curves for (K) and (K) in Figs. 9a,b show the spectral index because log10(-n) = log10((K)) - log10((K)). Our spectral index is nearly constant between the maxima in Figs. 9a,b, whereas it increases above and below the maxima. Our spectral index varies mostly with changes in , which is shown in Fig. 6.

Royrvik and Smith (1984) express reflectivity in terms of -n (K) such that (k1) from rocket data can be used to calculate reflectivity. Ulwick et al. (1988) used that formula to determine S/N of the Poker Flat radar for STATE; they use the rocket data for (k1) to calculate reflectivity and hence S/N, and compare it with the radar measurements of S/N. The portion of their figure 13a showing S/N in the turbulent layer in question is transcribed onto Fig. 10. Considering that their spectral level evaluated at the Bragg wave number has two maxima (see Fig. 9a or 9b), it is not known why only the lower maximum is seen in Fig. 10.

We performed a weighted running average of our dissipation rates to match the radar-range resolution used during STATE. That range resolution was approximately a Gaussian range profile having 300 m between full-width half-power points. We truncate the Gaussian at -250 m and +250 m from the center of a range-resolution volume; only 5% of the area under the Gaussian curve is beyond these truncation points. As for Figs. 5 and 6, the dissipation rates are also averaged over horizontal dimensions 1.13 by 1.13 km corresponding to the width of the radar beam at 88.75 km. We use the resultant dissipation rates in our turbulence advection model to calculate (K) at 100-m height intervals; next we calculate the cross section from Eq. (23), obtain reflectivity, 4, and use the relation by Ulwick et al. (1988) to calculate S/N for the Poker Flat radar. Our result is presented in Fig. 10 for comparison with the radar data and rocket-based calculation. Both upper and lower maxima of the simulation's curves in Figs. 9a,b are also evident in Fig. 10. In Fig. 10, the narrower height range within which the rocket-based S/N is large as compared to the simulation's S/N reflects this same feature in the spectral levels in Figs 9a,b. Agreement of rocket-based S/N with the simulation's S/N in Fig. 10 is a result of the fact that the ion models were purposefully chosen to fit the rocket-measured spectrum in Fig. 8. The lower edge of the layer in the radar's S/N appears about 400 m higher in height compared with the simulation's S/N layer, although both layers have comparable depth. Perhaps the shift in height is caused by the fact that the rocket trajectory was, for the height interval shown, about 9 km NNE of the center of the radar beam combined with the fact that we chose the simulation's height based on the rocket's profile in Fig. 4. The radar's wind data presented by Fritts et al. (1988) shows, at the height of the layer and the time of the rocket launch, a wind vector of 25 m s-1 toward the NNE. Thus, the separation of the rocket trajectory and radar beam is about 2.5 billow lengths in the streamwise direction. We expect streamwise coherence of KH billows over at least 3 billow lengths, so it is possible that the rocket pierced the turbulent layer about 1/2 billow length out of spatial phase from the center of the radar beam. On the other hand, at the wind speed of 25 m s-1, 1/2 billow length is advected in the streamwise direction in slightly more than a minute, and, for the radar's S/N shown in Fig. 10, the time at which that data was obtained and that data's averaging duration is not given by Ulwick et al. (1988) or Kelley and Ulwick (1988). The radar data does not show the two maxima that correspond to the top and bottom of the turbulent layer; the simulation's maxima would be reduced if a later simulation time was chosen. Clearly, we have uncertainties with regard to comparison of our simulation with the radar data that are in addition to the uncertainties mentioned previously with regard to comparison with the rocket data.

9. Discussion

Numerical solutions of the hydrodynamic equations are useful for investigation of turbulence processes in the upper atmosphere. Indeed, at altitudes near the mesopause, present computational resources suffice to realistically simulate KH dynamics throughout the full range of spatial scales from billow scale to viscous scale. We use these solutions to predict the morphology of turbulent layers and to calculate dissipation rates. The theory of multipolar diffusion is used to obtain the continuity equation of the ion mixing ratios. This equation is linearized and used to derive equations for spatial spectra and cospectra of ion mixing ratios by use of an average over a locally isotropic ensemble. A turbulence closure that agrees with data is applied to these equations. The dissipation rates from the hydrodynamic solution give the inertial-range spectral and cospectral levels that initialize the solution of the turbulence advection model toward high wave numbers. The electron mixing-ratio spectrum is then calculated from charge neutrality and the 3D electron number-density spectrum is calculated, from which the radar cross section is obtained.

The morphology of the STATE electron number-density fluctuations in Fig. 4 is similar to that of the potential temperature from the hydrodynamic solution despite the fact that electron number density is not a conserved quantity. The 1D electron spectrum is calculated and compared to STATE rocket measurements. Radar signal-to-noise ratio is calculated and compared to the radar measurements during STATE. Our results support the hypothesis that very massive ions mixed by turbulence is a reasonable mechanism for PMSE when turbulence is present.

Our present efforts are to further improve comparisons with data. We will calculate realistic radar scattering volumes from antenna theory, calculate reflectivity as a function of position within the scattering volume, and use that reflectivity field to calculate the first three Doppler spectral moments (power, velocity, width) by means of reflectivity weighted velocities. We will include gravity-wave breaking simulations, and consider initial ion mixing-ratio profiles that are more general than the linear profile, including variable ion composition in the vertical. Our future comparisons with data will include other experiments. Our calculations using the turbulence advection model will include effects of turbulence nonstationarity as well as the effects of coupling to potential-temperature fluctuations that arise, for instance, from the temperature dependence of the ion diffusivities.

It seems necessary to place our present theory in the context of the earlier work on turbulence advection applied to PMSE by Cho et al. (1992, 1996). We use the momentum, continuity, and Maxwell's equations (Eqs. 5-7) for each species of ion to derive our results. In contrast, Cho et al. (1992, 1996) assume that an effective electron diffusivity applies. In Sec. 5, we noted that an effective electron diffusivity does not apply. In fact (Hill, 1978a), the ionization has as many diffusion modes as there are distinct species of ions, and the behavior of the electron number density can differ greatly from that supposed on the basis of an effective electron diffusivity. Cho (1993) used computations of multipolar diffusion theory to assign values to an effective electron diffusivity for cases for which such a diffusivity is not valid; in fact, Cho (1993) observed more than one diffusion mode in his computations. One misunderstanding by Cho et al. (1992, 1996), which is caused by those assigned values, is that the Schmidt number can not be large unless the heavy-ion charge number density is a significant fraction of the electron number density. Actually, Schmidt numbers are determined from the momentum-transfer collision frequency; as such, Schmidt numbers are independent of charge number densities. Both we and Cho et al. (1992, 1996) use the turbulent advection model originated by Hill and Bowhill (1976) and subsequently refined.

Acknowledgments. We thank Kenneth A. Mitton and William D. Otto for their expert programming. This research was supported by the National Science Foundation under grant ATM-9618004 and the Air Force Research Laboratory under grant F19628-98-C-0030.

References

Andreassen, O., P.O. Hvidsten, D. C. Fritts, and S. Arendt, Vorticity dynamics in a breaking

gravity wave, 1. Initial instability evolution, J. Fluid Mech., 367, 27-46, 1998.

Arendt, S., D. C. Fritts, and O. Andreassen, The initial value problem for Kelvin vortex waves,

J. Fluid Mech., 344, 181-212, 1997.

Bjorn, L. G. , E. Kopp, U. Herrmann, P. Eberhardt, P. H. G. Dickinson, D. J. Mackinnon,

F. Arnold, G. Witt, A. Lundin, and D. B. Jenkins, Heavy ionospheric ions in the formation

process of noctilucent clouds, J. Geophys. Res. 90, D5, 7985, 1985.

Cho, J. Y. N., Radar scattering from the summer polar mesosphere: Theory and observations,

204pp., Ph.D. thesis, Cornell University, Ithaca, New York, 1993.

Cho, J. Y. N., T. M. Hall, and M. C. Kelley, On the role of charged aerosols in polar mesosphere

summer echoes, J. Geophys. Res., 97, 875-886, 1992.

Cho, J. Y. N., C. M. Alcala, M. C. Kelley, and W. E. Swartz, Further effects of charged aerosols

on summer mesosphere radar scatter, J. Atmos. Terr. Phys., 58, 661-672, 1996.

Cho, J. Y. N. and M. C. Kelley, Polar mesosphere summer radar echoes: observations and

current theories, Rev. Geophys, 31, 243-265, 1993.

Cho, J. Y. N. and J. Rottger, An updated review of polar mesosphere summer echoes:

observation, theory, and their relationship to noctilucent clouds and subvisible aerosols,

J. Geophys. Res., 102, 2001-2020, 1997.

Fritts, D. C., S. A. Smith, B. B. Balsley, C. R. Philbrick, Evidence of gravity wave saturation and

local turbulence production in the summer mesosphere and lower thermosphere during the

STATE experiment, J. Geophys. Res., 93, 7015-7025, 1988.

Fritts, D. C., S. Arendt, and O. Andreassen, Vorticity dynamics in a breaking internal gravity

wave, 2. Vortex interactions and transition to turbulence, J. Fluid Mech., 367, 47-65, 1998.

Fritts, D. C., T. L. Palmer, O. Andreassen, and I. Lie, Evolution and breakdown of Kelvin-

Helmholtz billows in stratified compressible flows, I: Comparison of two- and three-

dimensional flows, J. Atmos. Sci., 53, 3173-3191, 1996.

Havnes, O., J. Troim, T. Blix, W. Mortensen, L. I. Naesheim, E. Thrane, and T. Tonnesen, First

detection of charged dust particles in the Earth's mesosphere, J. Geophys. Res. 101,

10839-10847, 1996.

Havnes, O., L. I. Naesheim, T. W. Hartquist, G. E. Morfill, F. Melandso, B. Schleicher, J. Troim,

T. Blix, and E. Thrane, Meter-scale variations of the charge carried by mesospheric dust,

Planet Space Sci. 44, 1191-1194, 1996.

Hill, R. J., Nonneutral and quasi-neutral diffusion of weakly ionized multiconstituent plasma,

J. Geophys. Res., 83, 989-998, 1978a.

Hill, R. J., Models of the scalar spectrum for turbulent advection, J. Fluid Mech., 88, 541-562,

1978b.

Hill R. J. and S. A. Bowhill, Small-scale fluctuations in D-region ionization due to

hydrodynamic turbulence, Aeronomy Report No. 75, University of Illinois, Urbana, Illinois,

Nov. 1976.

Hill R. J. and S. A. Bowhill, Transient compressional response of D-region ionization, J. Atmos.

Terr. Phys., 39, 333-346, 1977.

Hill, R. J., and K. A. Mitton, Turbulence-induced ionization fluctuations in the lower ionosphere,

NOAA Technical Report ERL 454-ETL 68, November 1998. (Available from the author or

the National Technical Information Service, 5285 Port Royal Road, Springfield, VA, USA.)

Hocking, W. K., On the extraction of atmospheric turbulence parameters from radar backscatter

Doppler spectra-I. Theory, J. Atmos. Terr. Phys., 45, 89-102, 1983.

Hocking, W. K., An assessment of the capabilities and limitations of radars in measurements of

upper atmosphere turbulence, Adv. Space Res., 17, (11)37-(11)47, 1996.

Inhester, B., J. C. Ulwick, J. Cho, M. C. Kelley, and G. Schmidt, Consistency of rocket and radar

electron density observations: Implication about the anisotropy of mesospheric turbulence,

J. Atmos. Terr. Phys., 52, 855-873, 1990.

Kelley, M. C. and J. C. Ulwick, Large- and small-scale organization of electrons in the

high-latitude mesosphere: implications of the STATE data, J. Geophys. Res., 93, 7001-7008,

1988.

Klaassen, G. P., and W. R. Peltier, The onset of turbulence in finite-amplitude Kelvin-Helmholtz

billows, J. Fluid Mech. 227, 1-35, 1985.

Klostermeyer, J., On the formation of electron depletions at the summer polar mesosphere,

Geophys. Res. Lett. 23, 335-338, 1996.

Langevin, M. P., Une formule fondamentale de theorie cinetique, Annales de Chimie et de

Physique, series 8, 5, 245-288, 1905.

Lubken, F.-J., G. Lehmacher, T. Blix, U.-P. Hoppe, E. Thrane, J. Cho, and W. Swartz, First

in-situ observations of the Schmidt number within a PMSE layer, Geophys. Res. Lett. 20,

2311-2314, 1993.

Lubken, F.-J., K.-H. Fricke, and M. Langer, Noctilucent clouds and the thermal structure near the

Arctic mesopause in summer, J. Geophys. Res. 101, 9489-9508, 1996.

Lubken, F.-J., M. Rapp, T. Blix, E. Thrane, Microphysical and turbulent measurements of the

Schmidt number in the vicinity of polar mesosphere summer echoes, Geophys. Res. Lett. 25,

893-896, 1998.

Lumley, J. L., and H. A. Panofsky, The Structure of Atmospheric Turbulence, 239pp.,

Interscience Publishers, New York, 1964.

Muschinski, A., and C. Wode, First in situ evidence for coexisting submeter temperature and

humidity sheets in the lower free troposphere, J. Atmos. Sci., 55, 2893-2906, 1998.

McDaniel, E. W., Collision Phenomena in Ionized Gases, 775pp., Wiley Series in Plasma Physics, John Wiley & Sons, Inc., New York, 1964.

Ottersten, H., Radar backscattering from the turbulent clear atmosphere, Radio Sci. 4, 1251-1255,

1969.

Palmer, T. L., D. C. Fritts, and O. Andreassen, Evolution and breakdown of Kelvin-Helmholtz

billows in stratified compressible flows, II: Instability structure, evolution, and energetics,

J. Atmos. Sci., 53, 3192-3212, 1996.

Pedersen, A., J. Troim, and J. A. Kane, Rocket measurements showing removal of electrons

above the mesopause in summer at high latitude, Planet. Space Sci. 18, 945-947, 1970.

Rottger, J., Middle atmospheric studies with the EISCAT radars: polar mesosphere summer

echoes, pp. 369-387, Kluwer Academic Publishers, Netherlands, 1993.

Rottger, J., Polar mesosphere summer echoes: Dynamics and aeronomy of the mesosphere,

Adv. Space Res., 14, (9)123-(9)137, 1994.

Royrvik, O., and L. G. Smith, Comparison of mesospheric VHF radar echoes and rocket probe

electron number density measurements, J. Geophys. Res. 89, 9014-9022, 1984.

Tennekes, H., and J. L. Lumley, A First Course in Turbulence, 300 pp., MIT Press, Cambridge,

Mass., 1972.

Ulwick, J. C., K. D. Baker, M. C. Kelley, B. B. Balsley, and W. L. Ecklund, Comparison of

simultaneous MST radar and electron density probe measurements during STATE,

J. Geophys. Res., 93, 6989-7000, 1988.

Villars, F., and V. F. Weisskopf, On the scattering of radio waves by turbulent fluctuations of the

atmosphere, Proc. IRE 43, 1232-1238, 1955.

Watkins, B. J., C. R. Philbrick, and B. B. Balsley, Turbulence energy dissipation rates and inner

scale sizes from rocket and radar data, J. Geophys. Res., 93, 7009-7014, 1988.

Werne, J., and D. C. Fritts, Stratified shear turbulence: Evolution and statistics, Geophys. Res.

Lett., 26, 439-442, 1999.

______________________________________

R. J. Hill (e-mail: rhill@etl.noaa.gov), D. C. Fritts (e-mail: dave@colorado-research.com)



List of Figures



Fig. 1 Cross sections in a streamwise vertical plane displaying the 3D evolution of vorticity magnitude. The 8 nondimensional times span the dynamical evolution, including early quasi-2D billow roll-up (t = 49), transition to 3D structure consisting of streamwise-aligned convective rolls (t= 66 and 83), creation of a turbulent layer (t 83 to 164), and eventual restratification and stabilization of the flow (t > 200).



Fig. 2 The evolution of the potential-temperature field shown in the same cross section and at the same times as in Fig. 1.



Fig. 3 From right to left in Eq. (21) the 3 terms are the long-dashed, short-dashed, and dotted curves, respectively, and the solid curve is their sum, (k). The dotted curve lies above the solid curve because the middle term in Eq. (21) is negative. The diamond at the bottom of the graph marks the Bragg wave number, and from left to right the triangles mark kK, and the 3 values of kB for Sc = 3.24, 126, and 565, respectively. The case considered is that in Sec. 8 having a light positive ion and dissipation rates averaged over the height range 87.75-89.33 km at the edge of the KH billow.



Fig. 4 The dimensionalized potential temperature along the vertical line at the center of the computational domain (solid curve) and the electron number density (dashed curve) transcribed from figure 12 by Ulwick et al. (1988).

Fig. 5 Profiles of the dimensionalized dissipation rate of potential temperature from the KH simulation. The dissipation has been averaged over volumes that are 1.13 by 1.13 km in the horizontal and 100 m in the vertical directions. Open and filled circles correspond to averaging volumes positioned at the center and edge of the KH billow, respectively.



Fig. 6 Profiles of the dimensionalized energy-dissipation rate from the simulation. The dissipation has been averaged over volumes that are 1.13 by 1.13 km in the horizontal and 100m in the vertical directions. Open and filled circles correspond to averaging volumes positioned at the center and edge of the KH billow, respectively.



Fig. 7 1-D electron number-density power spectra from the turbulence advection model for 4 cases: edge and light positive ion (long-dashed curve); edge and heavy positive ion (dotted curve); center and light positive ion (short-dashed curve); center and heavy positive ion (solid curve).



Fig. 8 The electron number-density 1-D spectrum from the turbulence advection model (same as the long-dashed curve in Fig. 7) compared with that from the rocket measurements in the altitude range 88-89 km. The model spectrum is the long-dashed curve that becomes white dashes where it overlays the rocket data. The vertical dashed line indicates the Bragg wave number of the Poker Flat radar.





Fig. 9a For the heavy-positive-ion model, profiles of electron number-density 1-D spectra, (K) (diamonds), and 3-D spectra , (K) (circles), evaluated at the Bragg wave number K. Filled symbols, center of computational domain; open symbols, edge of domain. Dashed curve, STATE 1 data for (K).



Fig. 9b Same as Fig. 9a, except the light-positive-ion model is used.



Fig. 10 A portion of figure 13a by Ulwick et al. (1988) showing S/N from the radar data (long-dashed curve), and from their calculation of S/N based on the rocket-measured electron number-density spectrum (short-dashed curve). S/N from the simulation for 4 cases: edge and light positive ion (open diamonds); edge and heavy positive ion (open circles); center and light positive ion (filled diamonds); center and heavy positive ion (filled circles).

List of tables

Table 1. Parameters of the ionization for comparison of the turbulence advection model with STATE data. q/e is the ratio of the ion's charge to the magnitude of the electron's charge.

N /Ne (=Q /Qe) is the ratio of ion to electron number densities. The scale heights, H, are related to the gradients, G(z0), by H = -Q(z0) [G(z0)]-1, where z0 = 88.75 km. The number of water molecules attached is given as #H2O from which follows the radius, r (in nm), on the basis of the density of ice, and the Schmidt number, Sc, follows from Langevin's (1905) collision model.



q/e N/Ne H (km) #H2O r (nm) Sc


model with a heavy positive ion:

1 -1 1.5 1.0 30000 6.2 257
2 +1 2.5 1.96 15000 4.9 164


model with a light positive ion:

1 -1 1.5 1.0 10000 4.3 126
2 -1 0.5 1.0 100000 9.2 565
3 +1 3.0 1.69 15 0.49 3.24