Turbulent Breakdown of Internal

Atmospheric Waves: I The Theory of Atmospheric Waves

R. G. Stockwell

Department of Physics, The University of Western Ontario London, Canada N6A 5B7 Phone: 519-661-3929, Fax: 519-661-3129, Email: stockwell@uwo.ca

December 19, 1995 DRAFT

Abstract

This paper will examine the role that internal bouyancy waves, or gravity waves have on the atmosphere at large, and how these effects are realized. These gravity waves transfer energy and momentum from the troposphere to the middle atmosphere. The deposition of momentum and energy is through the turbulent breakdown of the gravity waves. Gravity waves have noticeable effects on the middle atmosphere in a number of ways:

• the pole-equator temperature difference is half what it should be (as predicted by purely radiative equilibrium).

• the maximum ozone density is at high latitudes rather than the tropics (which is predicted by photochemical equilibrium).

• the middle atmosphere (50 - 100 km in altitude) is warmer in winter than in summer due to the energy transported there (and deposited by turbulent breakdown) by gravity waves.

• the winds in the middle atmosphere blow ”the wrong way”, (ie opposite in direction to what a simple momentum budget would suggest) due to momentum dumped by breaking gravity waves.

This paper will survey the recent advances in understanding gravity waves, and the mechanisms in which they break down.

I. Introduction

Atmospheric wave motions are essential components of the general circulation of the atmosphere. Gravity waves (which are largely excited in the planetary boundary layer and propagate to altitudes upwards of 100km) have a profound effect on the large scale flow by producing turbulent dissipation that modifies the energy and momentum distributions.

Fig. 1. The temperature profile of the atmosphere [Hines, 1974]

Beacuse of the anisotropy due to gravity, the atmosphere has an exponential decrease in density with height (ie ρ = exp{−z/H}), with a scale height H of about 5 km. Due to the exponential decrease of density with height, the gravity waves grow in amplitude as they propagate upwards, until at some threshold the atmosphere can no longer support such large oscillations. At this critical layer the wave is said to break. In doing so, it deposits energy into that atmospheric region which manifests itself in an increase in turbulence.

In short, these internal waves are created in turbulent conditions (source regime), propagate a large distance away without permanently affecting the medium through which it is travelling (linear wave regime), and then deposit their energy and momentum when the wave breaks (absorption regime).

II. Gravity Waves

A theoretical treatment of atmospheric disturbances was succesfully made by Hines (1960). It is based on an idealized model of the atmosphere where the atmosphere is isothermal, windless, and inviscid, and where the waves themselves are small perturbations on the background (ie. only low order terms are kept). Since the characteristic sizes of the waves are small (and the Coriolis force is small in comparison to the force of gravity), the Coriolis force is neglected in the following derivation. That this is a valid approximation is borne out in a comparison of the theoretical predictions and experimental observations.

The main equations are:

The equation of motion where the forces considered are the pressure gradient term and the gravitational force.

\frac{DU}{Dt} = −\frac{1}{ρ}∇p + g  (1)

where \frac{DU}{Dt} is the Stoke’s or advective derivative.

\frac{DU}{Dt} = \frac{∂U}{∂t} + \frac{1}{2}∇U^{2} + (∇ × U) × U  (2)

As is common in the literature, this expression is written in a shorthand form (in analogy with the advective derivitive of a scalar) so that the equation of motion becomes

\frac{∂U}{∂t} + (U · ∇)U = \frac{−1}{ρ}∇p + g  (3)

The equation of mass continuity

\frac{∂ρ}{∂t} + ∇ · (ρU) = 0  (4)

The adiabatic equation which is also the equation for energy conservation

\frac{∂p}{∂t} + U · ∇p = c^{2}\{\frac{∂ρ}{∂t} + U · ∇ρ\}  (5)

This model of the atmosphere allows for gravity waves as well as acoustic waves where the forces involved are due to compressibility of the atmosphere. By symmetry in the horizontal plane, without loss of generality one can take the horizontal components to be in one direction along the x-axis.

One then assumes that the basic quantites (p, ρ, and U) can be written as zero and first order terms:

p = p_{o} + p_{1}

ρ = ρ_{o} + ρ_{1}
U = U_{o} + U_{1}  (6)

Substituting the above into the equations of motions, continuity and state, we then get the zero order equations which represent the variations in the background quantities:

\frac{DU_{ox}}{Dt} = −\frac{1}{ρ_{o}} \frac{∂p_{o}}{∂x}  (7)

which is the gradient wind equation in the absence of Coriolis forces,

\frac{DU_{oz}}{Dt} = −\frac{1}{ρ_{o}} \frac{∂p_{o}}{∂z} + g  (8)

and providing that Uo = 0 (ie, the no-wind case)

\frac{∂p_{o}}{∂t} = c^{2} \frac{∂ρ_{o}}{∂t}  (9)

and

\frac{∂ρ_{o}}{∂t} = 0  (10)

The last equation tells us that the background density may be treated as time independant. It is however, still a very strong function of altitude z.

The conditions for the isothermal, windless, flat atmosphere are

\frac{∂T}{∂z} = \frac{∂T}{∂x} = U_{ox} = U_{oz} = 0  (11)

Substituting this into the first order equations of motions, continuity and state, we get

ρ_{o} \frac{∂U_{1x}}{∂t} = −\frac{∂p_{1}}{∂x}  (12)
ρ_{o} \frac{∂U_{1z}}{∂t} = −\frac{∂p_{1}}{∂z} − gρ_{1}  (13)
\frac{∂p_{1}}{∂t} + U_{1z} \frac{∂p_{o}}{∂z} = c^{2}\{\frac{ρ_{1}}{∂t} + U_{1z} \frac{∂ρ_{o}}{∂z} \}  (14)
\frac{∂ρ_{1}}{∂t} + U_{1z} \frac{∂ρ_{o}}{∂z} + ρ_{o}\{\frac{U_{1x}}{∂x} + \frac{∂U_{1z}}{∂z} \} = 0  (15)

This set of equations governs the behaviour of linear waves in a flat isothermal atmosphere. When one substitutes the following wavelike form into the governing equations

\frac{p_{1}}{p_{o}P} = \frac{ρ_{1}}{p_{o}R} = \frac{U_{1x}}{X} = \frac{U_{1y}}{Y} = \frac{U_{1z}}{Z}

= A_{o}exp{i(ωt − K_{x}x − K_{y}y − K_{z}z)}  (16)

one gets the following dispersion relation

ω^{4} − ω^{2}c^{2}(K_{x}^{2} + K_{z}^{2}) + (γ − 1)g^{2}K_{x}^{2} + iω^{2}γgK_{z} = 0  (17)

where we will use the convention that complex values are capitalized and real values are lower case, ie. wavenumber K = k + i · {Im[K]}. By solving the quadratic equation in ω^{2} one sees that both K_{x} and K_{z} cannot be real and non-zero. Thus the wave must have attenuation or growth in either the vertical or the horizontal direction. By symmetry in the horizontal plane, one can assume that there is no variation in amplitude in the horizontal plane. Thus K_{x} = k_{x} and it is purely real.

The real part of the dispersion relation becomes

ω^{4} − ω^{2}c^{2}(k_{x}^{2} + k_{z}^{2} − (Im[K_{z}])^{2}) − ω^{2}γgIm[K_{z}] + (γ − 1)g^{2}k_{x}^{2} = 0  (18)

and the imaginary part becomes

ω^{2}γgRe[K_{z}] − 2ω^{2}c^{2}Re[K_{z}]Im[K_{z}] = 0  (19)

so providing that the frequency ω is not zero, either

Re[K_{z}] = 0  (20)

or

Im[K_{z}] = \frac{γg}{2c^{2}} = \frac{1}{2H}  (21)

(where H is the scale height H = RT/Mg and R is the gas constant). The first condition means that the vertical wavenumber is purely imaginary, and we have attenuation or growth of the wave’s amplitude in the vertical direction (which is a characteristic of evanescent and surface waves). This category of waves is referred to as external waves. It is apparent from the above equations that external waves exhibit no vertical phase variation, while the internal waves do. Thus the origin of the term internal atmospheric gravity waves.

There are four solutions to the fourth order equation, comprising two distinct regimes. One regime is where the frequency is greater than ω_{a} defined by

ω_{a} = \frac{γg}{2c} = \frac{c}{2H}  (22)

This frequency ω_{a} is known as the acoustic frequency. The waves with frequencies larger than the acoustic frequency are called acoustic waves.

The other regime consists of the lower frequency waves, which has an upper bound ω_{g} called the bouyancy frequency defined by

ω_{g}^{2} = \frac{(γ − 1)g^{2}}{c^{2}}  (23)

This is shown below to be the frequency of oscillation of an air parcel displaced in a stably stratified atmosphere known as the Brunt-Vaisala frequency N. These long-period waves are referred to as gravity waves. The four solutions to the dispersion relation refer to two waves of each type with opposite direction of phase propagation.

Fig. 2. The possible modes for gravity waves in the atmosphere as a function of vertical wavenumber k_{z} [Beer, 1974]

Fig. 3. The possible modes for gravity waves in the atmosphere as a function of horizontal wavenumber k_{h} [Beer, 1974]

If the Coriolis term had been kept, yet another solution would be possible. These would correspond to Rossby or planetary waves, with wavelengths of thousands of kilometers, and which drive the local weather patterns at mid-latitudes. These waves, however, are beyond the scope of this paper.

One can rewrite the dispersion relation in terms of the acoustic ω_{a} and bouyancy ω_{g} frequencies

k_{z}^{2} = k_{x}^{2}(\frac{ω_{g}^{2}}{ω^{2}} − 1) + \frac{(ω^{2} − ω_{a}^{2})}{c^{2}}  (24)

or

k_{x}^{2} = \frac{ω^{2}}{c^{2}} (\frac{ω^{2} − ω_{a}^{2} − c^{2}k_{z}^{2}}{ω^{2} − ω_{g}^{2}} )  (25)

We can then rearange them into the form

\frac{k_{x}^{2}}{l^{2}} + \frac{k_{z}^{2}}{m^{2}} = 1  (26)

where

l^{2} = (\frac{ω^{2} − ω_{a}^{2}}{ω^{2} − ω_{g}^{2}} )\frac{ω^{2}}{c^{2}}  (27)

and

m^{2} = (\frac{ω^{2} − ω_{a}^{2}}{c^{2}} )  (28)

to put them into the form of a general conic. In the real atmosphere ω_{a} is always greater than ω_{g} since γ < 2. For acoustic waves l^{2} > 0 and m^{2} > 0, while for gravity waves l^{2} > 0 and m^{2} < 0. For the special case where ω_{a} > ω > ω_{g} the dispersion relation is

\frac{k_{x}^{2}}{|l|^{2}} + \frac{k_{z}^{2}}{|m|^{2}} = −1  (29)

and we have no real valued solution.

It is illuminating to rewrite this dispersion relation in terms of the two limiting frequencies:

ω^{2} = f^{2}sin(θ) + N^{2}cos(θ)  (30)

where θ = atan(\frac{k_{z}}{k_{x}}), (ie. the angle of propagation of the wave) and N is the Brunt-Vaisala frequency and f is the Coriolis parameter.

A. The Coriolis Parameter and the Brunt-Vaisala Frequency

As can be seen above, linear gravity waves can exist between two fundamental frequencies in the atmosphere. On the low frequency end the particle displacement is horizontal and the Coriolis force is the dominant term (Figure 4a).

8 TURBULENCE AM551 TERM PAPER

Fig. 4. The particle displacement for the two limiting frequencies f and N [Garret, 1970]

The Coriolis parameter is derived from the Coriolis force on a mass m travelling with velocity U

|2m(ΩE × U)| = 2mΩ_{E}|U|cos(θ)

= mf |U|  (31)

where θ is the latitude. Thus,

f = 2Ω_{E}cos(θ)  (32)

where Ω_{E} = 7.29 × 10^{−5} s^{−1} is the rotation frequency of the earth.

Fig. 5. The components of rotation at a point on the earth’s surface

For mid-latitudes f corresponds to a period of about 13 hours.

DRAFT December 19, 1995

The Brunt-Vaisala frequency is the natural frequency of oscillation of an air parcel and is a measure of stability of the atmosphere. If a parcel of air is displaced vertically in a stable atmosphere, it will feel a restoring force (Figure 6). If the displacement is small, then the forces due to gravity are comparable to those due to compression. If the displacement was large, and compression forces dominated over the gravity force and the wave would propagate as an longitudinal acoustic wave.

Fig. 6. The bouyancy of an air parcel in the atmosphere [Sica, 1995]

Consider a parcel of air at height z. If it is displaced vertically a distance Δz, then the equation of motion is

ρ_{o} \frac{d^{2}Δz}{dt^{2}} = g(Δρ_{o} − Δρ)  (33)

For the first term on the right hand side, we can write

Δρ_{o} = (\frac{dρ_{o}}{dz})Δz  (34)

For the second term on the right hand side, we can use the hydrostatic approximation

\frac{Δp}{Δz} = −ρ_{o}g  (35)

and the fact that for an adiabatic process the speed of sound is given by

c^{2} = \frac{Δp}{Δρ}  (36)

to rewrite the equation of motion as

ρ_{o} \frac{d^{2}Δz}{dt^{2}} = g(\frac{dρ_{o}}{dz} + \frac{ρ_{o}g}{c^{2}} )Δz  (37)

This is in the form for a simple harmonic oscillator and we can see that provided

\frac{g}{ρ_{o}} (\frac{dρ_{o}}{dz} + \frac{ρ_{o}g}{c^{2}} ) < 0  (38)

the parcel will oscillate with an angular frequency

N^{2} = −g[\frac{d(lnρ)}{dz} + \frac{g}{c^{2}} ]  (39)

The term in equation [38] is a measure of static stability. If it was greater than zero than we would have had an unstable equilibrium, and a displaced parcel of air would continue to move away from its starting position. (Likewise, for a downward displacement the air parcel would continue to accelerate downwards as the gravitational force would be stronger than the hydrostatic force.) Typical values for the Brunt-Vaisala frequency correspond to a period of about 5 to 10 minutes.

III. Sources of Gravity Waves

Gravity waves can be generated by a wide range of forcing events. Many (but by no means all) are generated in the tropophere and propagate upwards. Because of the increase in amplitude of the waves with height, it is these waves that are primarily of interest.

The primary sources of gravity waves include the following:

1. Topography (orographic forcing) The problem of air flow over a mountain has been studied for many years by people interested in the boundary layer effects and turbulence, as well as people interested in the generation of gravity waves. The perturbation to the flow caused by the mountain gives rise to motions with scales and initial amplitudes determined by the terrain itself.

Fig. 7. Generation of gravity waves by orographic forcing [Beer, 1974]

An excellent example of gravity waves caused by orographic forcing are the mountain waves ocurring on the lee side of mountains (Figure 7).

2. Convection Convective instabilites can cause convective plumes to penetrate upwards into stably stratified regions [Ahrens, 1988]. The plumes expand and cause the air layer to oscillate causing gravity waves. A weather front is a primary location for large scale convection, and a good source for gravity waves.

3. Polar front jet stream The strong velocity shear associated with the jet stream gives rise to dynamical instabilites which are a large source for gravity waves.

4. Breaking tidal and planetary waves As these larger scales structures undergo absorption, they degenerate into smaller scale waves. This is evidence of the energy cascade to smaller scales that is very common in atmospheric fluid dynamics.

5. Explosions Explosions, whether manmade or natural (volcanic eruptions) are a point source generator of gravity waves (Figure 8).

Fig. 8. The propagation of gravity waves from a point source (explosion) [Hines, 1960]

Above ground nuclear tests have been seen to produce gravity waves that propagated throughout the whole atmosphere. The gravity waves generated by a particular detonation (Figure 8) had periods that increased with distance from the source as predicted by Hines theory. If the gravity waves propagate up to the thermosphere, they can propagate for extremely large distances, and are called travelling ionospheric disturbances. These TIDs have been used to monitor above ground nuclear tests.

IV. Instability Mechanisms

Instability will refer to any situation in which energy is transferred from the wave to the mean flow. As is well known, one property of a linear wave is that is does not permanently displace the medium through which it propagates. If this were the case in the atmosphere, then the study of internal waves would interest only the curious few. But the amplitude of the wave (ie. the perturbations in wind due to the wave) grows exponentially with height. This growth with height of the wave can seen as a consequence of conservation of energy.

Fig. 9. The growth of a gravity wave’s amplitude with height [Hines,1960]

The kinetic energy density is given by

KE = \frac{1}{2}ρv^{2} = constant  (40)

therefore the square of velocity is inversely proportionally to density

v = \sqrt{\frac{2(constant)}{ρ}}  (41)

Since the density decreases as

ρ(z) = ρ_{0}exp(\frac{−z}{H})  (42)

where H is the scale height. Due to the exponential growth of the wave’s amplitude with height, however, instablities will form at some critical layer. At such a level, the linear wave becomes unstable and breaks. It is through these instabilities that the waves play a significant role in the determining the circulation of the middle atmosphere. Perhaps the most dramatic effects of gravity wave motions are the formation of the zonal jet streams, as well as the thermal structure of the upper atmosphere.

To start with an examination of a critical layer, we have to generalize the treatment leading to our governing equations for internal gravity waves. (It should be noted that evanescent waves have no vertical group or phase velocity since they do not propagate upwards.) First we have to introduce the case for a constant horizontal wind. It can be shown [Beer, 1974], that the result of applying a constant horizontal wind |U| = U_{ox} would result in a doppler shift of the frequency

Ω = ω − k_{x}U_{ox}  (43)

as one would expect.

If we also relax the isothermal atmosphere requirement to that of a slowly varying atmosphere, then we get the following relations [Beer, 1974]

ω → Ω = ω − Uo · k  (44)
ω_{g} → ω_{b} = \sqrt{ω_{g}^{2} + \frac{g}{T} \frac{∂T}{∂z}}  (45)
ω_{a} → ω_{an} = \sqrt{ω_{a}^{2} + \frac{γ}{2} \frac{g}{T} \frac{∂T}{∂z}}  (46)

Now we can substitute these into the dispersion relation for gravity waves (24)

k_{z}^{2} = k_{x}^{2}(\frac{ω_{g}^{2}}{ω^{2}} − 1) + \frac{(ω^{2} − ω_{a}^{2})}{c^{2}}  (47)

and differentiate it to give us the relation for group velocity in the vertical direction

V_{gz} = \frac{∂Ω}{∂k_{z}}

V_{gz} = \frac{Ω^{3}k_{z}c^{2}}{Ω^{4} − k_{x}^{2}N^{2}c^{2}}  (48)

As can clearly be seen from the above equations, as the wave’s horizontal phase velocity approaches the horizontal wind, ie

V_{px} = \frac{ω}{k_{x}} → U_{ox}  (49)

that the group velocity becomes zero. This occurs at the so called critical layer. It has been observed that very little wave energy propagates through a critical layer. Hines (1967) suggests that the condition Ω = 0 represents an absorption of the wave. Also note that as Ω → 0 that the vertical wavenumber k_{z} → ∞, so that the vertical scales size of the wave becomes small, and viscous dissipation increases.

Several mechanisms contribute to the breakdown of gravity waves. Along with convective and shear instabilities, there are other mechanisms such as molecular diffusion, radiative damping, inertial instabilites and the cascade of energy to small scales by nonlinear wave-wave interactions. This latter list operates most effectively on waves with small spatial scales and low group velocities. Thus spatial and temporal scales suggest that these processes are relatively less important than the shear and convective instabilities in the dissipation of the dominant wave motions seen in the atmosphere.

The following sections highlight the instability mechanisms that are frequently encountered in atmospheric studies.

A. Kelvin Helmholtz Instability

Studies of shear instabilities date back to the work of Helmholtz [1868], Kelvin[1871] and Rayleigh [1892]. Taylor [1931] studied heterogenous shear flows with uniform shear \overline{\frac{∂u}{∂z}} and static stability N^{2} = −\frac{g}{\overline{ρ}} \frac{∂\overline{ρ}}{∂z} , and found that the flow was stable for a Richardson number

Ri ≡ N^{2}/\overline{\frac{∂u}{∂z}}  (50)

The unstable wave motion examined in all the early papers was vertically evanescent above and below the shear layer, and is known as the Kelvin Helmholtz Instability (KHI). This has long been recognized as a significant mechanism for producing turbulence and mixing in stabily stratified fluids. KHI is one of the main mechanisms for gravity wave breakdown (along with convective instabilites to be described later).

The initial stages of the KH instability are well understood. Numerical computations by Patnaik (1976) and Peltier (1978) have confirmed linear stability predictions of the growth rates and structure of the growing waves. To gain a better insight to the mechanism of KH instability, the results of a Reynolds-stress turbulence second-order closure model on breaking KH billows will be presented. This section largely follows Sykes and Lewellen [1982].

One problem with the previously mentioned models is that the laminar calculations at moderate Reynolds numbers have been unable to identify the secondary instability which results in the breakdown of the vortex into a turbulent layer. By using a second order closure, Sykes is extending the analysis to study the breakdown of the rolled-up vortex layer and the subsequent spread of turbulence. This approach is thought to be of some validity because the model calculates the growth of turbulent correlations moderately accurately. It is expected to be more valid than a ’mixing length’ closure model since the time scale of the growth of turbulence is the same as the mean flow time scale, and the turbulence may not be in local equilibrium, as the mixing length model assumes.

The closure scheme is described in detail by Lewellen (1977) so we present only the final equations. Assuming an incompressible, Boussinesq flow in two dimensions (x, z) with z vertical, and defining a stream function ψ, such that U = \frac{∂ψ}{∂z} and W = −\frac{∂ψ}{∂x} then the vorticity is ζ = ∇^{2}ψ, and the equations of motion are:

\frac{∂ζ}{∂t} + U \frac{∂ζ}{∂x} + W \frac{∂ζ}{∂z} = −\frac{g}{T_{o}} \frac{∂T}{∂x} − \frac{∂^{2}}{∂z^{2}} \overline{uw} − \frac{∂^{2}}{∂x∂z} (\overline{uu} − \overline{ww}) + \frac{∂^{2}}{∂x^{2}} \overline{uw},  (51)
\frac{∂T}{∂t} + U \frac{∂T}{∂x} + W \frac{∂T}{∂z} = −\frac{∂}{∂x}\overline{uθ} − \frac{∂}{∂z} \overline{wθ},  (52)

\frac{D}{Dt}\overline{u^{2}} = −2\overline{uu}\frac{∂U}{∂x} − 2\overline{uw} \frac{∂U}{∂z} + v_{c} \frac{∂}{∂x}(qΛ\frac{∂\overline{u^{2}}}{∂x}) + v_{c} \frac{∂}{∂z} (qΛ\frac{∂\overline{u^{2}}}{∂z}) − \frac{q}{Λ}(\overline{uu} − \frac{q^{2}}{3} ) − \frac{2bq^{3}}{3Λ} ,  (53)
\frac{D}{Dt}\overline{v^{2}} = v_{c} \frac{∂}{∂x}(qΛ\frac{∂\overline{v^{2}}}{∂x}) + v_{c} \frac{∂}{∂z} (qΛ\frac{∂\overline{v^{2}}}{∂z}) − \frac{q}{Λ}(\overline{vv} − \frac{q^{2}}{3} ) − \frac{2bq^{3}}{3Λ} ,  (54)

\frac{D}{Dt}\overline{w^{2}} = −2\overline{uw} \frac{∂W}{∂x} − 2\overline{ww} \frac{∂W}{∂z} + \frac{2g}{T_{o}} \overline{wθ} + v_{c} \frac{∂}{∂x}(qΛ\frac{∂\overline{w^{2}}}{∂x}) + v_{c} \frac{∂}{∂z} (qΛ\frac{∂\overline{w^{2}}}{∂z}) − \frac{q}{Λ}(\overline{ww} − \frac{q^{2}}{3} ) − \frac{2bq^{3}}{3Λ} , (55)

\frac{D}{Dt}\overline{uw} = −\overline{uu}\frac{∂W}{∂x} − \overline{ww} \frac{∂U}{∂z} + \frac{g}{T_{o}} \overline{uθ} + v_{c} \frac{∂}{∂x}(qΛ\frac{∂\overline{uw}}{∂x}) + v_{c} \frac{∂}{∂z} (qΛ\frac{∂\overline{uw}}{∂z}) − \frac{q}{Λ}(\overline{uw}),  (56)

\frac{D}{Dt}\overline{uθ} = −\overline{uu}\frac{∂T}{∂x} − \overline{uw} \frac{∂T}{∂z} − \overline{uθ} \frac{∂U}{∂x} − \overline{wθ} \frac{∂U}{∂z} + v_{c} \frac{∂}{∂x}(qΛ\frac{∂\overline{uθ}}{∂x}) + v_{c} \frac{∂}{∂z} (qΛ\frac{∂\overline{uθ}}{∂z}) − \frac{Aq}{Λ} (\overline{uθ}),  (57)

\frac{D}{Dt}\overline{wθ} = −\overline{uw} \frac{∂T}{∂x} − \overline{ww} \frac{∂T}{∂z} − \overline{uθ} \frac{∂W}{∂x} − \overline{wθ} \frac{∂W}{∂z} + v_{c} \frac{∂}{∂x}(qΛ\frac{∂\overline{wθ}}{∂x}) + v_{c} \frac{∂}{∂z} (qΛ\frac{∂\overline{wθ}}{∂z}) − \frac{Aq}{Λ} \overline{wθ} + \frac{g}{T_{o}} \overline{θ^{2}}, (58)

\frac{D}{Dt}\overline{θ^{2}} = −2\overline{uθ} \frac{∂T}{∂x} − 2\overline{wθ} \frac{∂T}{∂z} + v_{c} \frac{∂}{∂x}(qΛ\frac{∂\overline{θ^{2}}}{∂x}) + v_{c} \frac{∂}{∂z} (qΛ\frac{∂\overline{θ^{2}}}{∂z}) − \frac{2bsq}{Λ} \overline{θ^{2}},  (59)

\frac{D}{Dt}Λ = 0.35\frac{Λ}{q^{2}} × (\overline{uu}\frac{∂U}{∂x} + \overline{uw} \frac{∂U}{∂z} + \overline{uw} \frac{∂W}{∂x} + \overline{ww} \frac{∂W}{∂z} )

+ 0.6bq + v_{c} \frac{∂}{∂x}(qΛ\frac{∂Λ}{∂x}) + v_{c} \frac{∂}{∂z} (qΛ\frac{∂Λ}{∂z}) − \frac{0.375}{q}[(\frac{∂qΛ}{∂x})^{2} + (\frac{∂qΛ}{∂z})^{2}] + \frac{0.8Λ}{q^{2}} \frac{g}{T_{o}} \overline{wθ}.  (60)

In the above, (U, W) is the mean velocity, T the mean temperature, (u, v, w) the mean turbulent velocity fluctuation, θ the temperature fluctuation, and q = [u^{2} +v^{2} +w^{2}]^{\frac{1}{2}} , Λ is a length scale associated with the energy containing eddies. The values of empirically determined constants are b = \frac{1}{8}, A = 0.75, v_{c} = 0.3, and s = 1.8.

For this calculation, Sykes started with an initial temperature profile given by

T = T_{o} + ΔT tanh\frac{z}{δ} ,  (61)

and an initial velocity profile given by

U = ΔUtanh\frac{z}{δ} .  (62)

These initial conditions define a Richardson number

Ri = (\frac{g}{T_{o}})(\frac{ΔT σ}{ΔU^{2}} )  (63)

and the calculation implicity has high a Reynold’s number because laminar diffusive terms have been excluded from the governing equations. Also, initial values for turbulence energy q^{2} and the turbulent length Λ scale were given. All other parameters were set to zero. periodic boundary conditions were assumed in the x direction and the boundaries in the z direction were placed far away from the shear layer so that they had no effect on the development of the instability.

The equations were solved by Sykes using “centered spatial derivitives, and a ADI scheme for the temporal derivatives” and the Poisson equation for the streamline function was solved directly “using the decomposition method of Swarztrauber and Sweet [1979]”.

For a Richardson number of 0.1, they initiated the KH instability by applying a very small sinusoidal variation in vorticity. The turbulence energy was set to q^{2} = 10^{−4}ΔU^{2} throughout the domain, and Λ was set equal to 0.4δ. In the following τ is a time parameter equal to τ = t(gΔT)/T_{o}ΔU).

Fig. 10. The contours for dimensionless temperature (T − T_{o})/ΔT h for the case Ri = 0.1. (a) τ = 3.0, (b) τ = 4.9, (c) τ = 5.9, (d) τ = 6.9, (e) τ = 8.9, (f) τ = 14.2 [Sykes, 1982]

With these parameters set, the billow developed as a laminar flow with a growing vortex core in the center, and thin braids forming between the cores. The shear is highest in these thin braids, and the Richardson number is lowest.

At τ = 4.9 the instability has almost reached its maximum amplitude. The core has become statically unstable and the q^{2} contours (Figure 11) show growth of turbulence associated with the convective instability.

At the next time step (Figure c) the temperature structure has been mixed by the turbulence which has reached a high level at the edge of the core. The core has begun to spread in the x direction.

Fig. 11. The contours for dimensionless turbulence kinetic energy q^{2}/ΔU^{2} for the same times as in Figure 10 [Sykes, 1982]

At τ = 8.9, the braids have mixed and turbulence has spread horizontally through the layer.

It is interesting to note that the convective instability is needed to create the turbulent breakdown of the KH instability, thus tying the two phenomenom closely together. This was a major feature missing from earlier studies of KHI. Also, small scale structure is seen in the temperature field, which is described in the paper [Sykes, 1982] as travelling wave-like disturbances moving along the top and bottom of the layer, thus implying that KHI are sources for gravity wave generation.

In his paper Sykes compare the results of his calculations with a tank experiment done by Thorpe [1973]. There was good qualitative agreement with the Thorpe results.

An examination of the energetics of the KHI is helpful to gain some insight. There is an initial buildup of roll energy as the growing billow extracts kinetic energy from the mean shear, while the turbulence energy changes very little. At about τ = 4.5, the turbulence energy begins to increase as the billow starts to break, and the kinetic energy hits its maximum. The kinetic energy is effectively transferred to turbulence energy via potential energy by the KHI.

Fig. 12. The cycle for the KH instability

The final value for the Richardson number that this model gets is between 0.35 and 0.40. The mixing process apparently extracts energy from the initial flow until it is unable to extract any more energy.

Fig. 13. An observation of a KH billow from August 6 1960 [Fritts 1985]

B. Convective Instability

When a parcel of air initially at rest and in equilibrium with the atmosphere absorbs heat it will expand. It becomes less dense than the surrounding cooler air and is bouyed upwards. This vertical exchange is called convection.

The rate of change of temperature with height that an (dry) air parcel undergoes when it is allowed to move adiabatically (without exchanging heat with the surroundings) is called the dry adiabatic lapse rate

Γ_{d} = −\frac{dT}{dZ} = \frac{g}{c_{p}}  (64)

where g is the acceleration due to gravity and c_{p} is the specific heat capacity of the atmosphere at constant pressure. This parameter is an indication of the convective stabilty of the atmosphere. If the vertical temperature gradient of the atmosphere −\frac{dT}{dZ} exceeds Γ_{d} then convection will take place and the layer becomes unstable. The convection would continue until the temperature gradient reaches the dry adiabatic lapse rate and becomes stable again.

V. Early Models of Atmospheric Circulation

Early attempts at modelling the circulation of the middle atmosphere ran into severe difficulties. It was found that small scale dynamics had a profound influence on the large scale circulation, and could not be ignored. In order to compensate for the discrepancies of the early models, a term called Rayleigh drag was introduced into the models that decelerated the mean flow so that the model compared more favourably with observations. In order to justify this, linear saturation theory was invented. In this theory a gravity wave grows in amplitude with height until it reaches a level in which the atmosphere can no longer sustain such growth. At this point the wave saturates and propagates with constant amplitude (Figure 14). The wave is now losing energy as it propagates higher, and this energy is lost in the production of turbulent eddies. It is these eddies that are said to be responsible for eddy diffusion and the Rayleigh drag.

This theory is incomplete because it has been found that the gravity wave field also accelerates the mean flow as well as applying a drag, and the interaction of the waves with the mean flow is much more complicated than the linear saturation just described.

VI. Spectrum of Middle Atmosphere Fluctuations

Several observation have been made on horizontal wavnumber power spectral densities which all have a common ’universality’. Originally it was proposed that the spectra was due to inertial range turbulence since it was very similar to the well known spectral form

Ψ(k_{x}) = c_{1}ϵ^{2/3}k_{x}^{−5/3}  (65)

An alternative explanation was that the fluctuations were due to gravity waves. The dynamical theory for the spectrum of internal waves has been studied for many years. Garret and Munk [1975] started this

Fig. 14. The linear saturation of an upwardly propagating gravity wave

effort with their theory of internal waves in the ocean. Their model showed that observations of ocean spectra would have a slope of -2.

Van Zandt [1982] suggested that this oceanic spectrum had a analogous form in the atmosphere. He showed that the power spectra of mesoscale wind fluctuations can be described by a universal spectrum due to internal atmospheric waves. The model of Garret and Munk for the ocean spectra could be applied to atmospheric horizontal velocity spectra with only minor modifications. These modifications include a change in the slope of the spectra from -2 to -5/3.

A. The Similitude Argument of Dewan

This section breifly outlines the argument of Dewan for the appearance of the universal spectrum in the atmosphere [Dewan, 1991].

Starting at Newton’s Second Law

ρ_{o} \frac{∂u}{∂t} = −\frac{∂P′}{∂x}  (66)
ρ_{o} \frac{∂v}{∂t} = −\frac{∂P′}{∂y}  (67)
ρ_{o} \frac{∂w}{∂t} = −\frac{∂P′}{∂z} − ρ′g  (68)

and the equation for incompressible fluid

\frac{∂u}{∂x} + \frac{∂v}{∂y} + \frac{∂w}{∂z} = 0  (69)

and

\frac{∂ρ′}{∂t} + w \frac{∂ρ_{o}}{∂z} = 0  (70)

Fig. 15. The horizontal wavenumber spectrum of gravity waves in the atmosphere [Sica, 1995]

By eliminating the density and pressure terms one gets [Gill, 1982]

\frac{∂^{2}(∇^{2}w)}{∂t^{2}} + N^{2}∇_{H}^{2}w = 0  (71)

where

N^{2} = −\frac{g}{ρ_{o}} \frac{dρ_{o}}{dz}  (72)

and

∇_{H}^{2} = \frac{∂^{2}}{∂x^{2}} + \frac{∂^{2}}{∂y^{2}}  (73)

When one makes the assumptions that we have wavelike solutions and that we have the usual dispersion relation for gravity waves, then we obtain the following conditions:

\frac{∂v}{∂z} = ik_{z}v  (74)

and

\frac{∂u}{∂z} = ik_{z}u  (75)

With these conditions and following linear inviscid theory where the necessary condition for dynamic instability is

Ri < \frac{1}{4}  (76)

then we can write [Phillips, 1977] that the condition on shear at the border line of turbulence is

S_{cr}^{2} = \{[\frac{∂v}{∂z}]^{2} + [\frac{∂u}{∂z}]^{2}\}_{max} = 4N^{2}  (77)

Choosing a coordinate system so that the horizontal axis is in the direction of horizontal velocity we can write this as

S_{cr}^{2} ≡ (\frac{∂u}{∂z})_{max}^{2} = 4N^{2}  (78)

Now, in order to examine the physical meaning of S_{cr} let

u = u_{o}sin(k_{z}z)  (79)

At critical shear then, we have

|(\frac{∂u}{∂z})_{max}| = |k_{z}u_{o}cosk_{z}z| = 2N  (80)

where the argument of the cosine must be 0^{o} or 180^{o} by the condition of being a maximum. For typical values of atmospheric parameters N = 2 × 10^{−2}s^{−1} and λ_{z} = 500m we get for the critical velocity

(u_{o})_{critical} = \frac{λ_{z}}{2π} (2N) = 3.18m/s  (81)

For velocities above this value, the instability will grow with time. There are several points to note. The thickness of the turbulent patch does not depend on λ. Also, the turbulent patch can be much smaller than the wavelength of the wave causing it.

So we have

(u_{o})_{cr}^{2} = \frac{4N^{2}}{k_{z}^{2}}  (82)

Using this and the fundamental fact that the total variance of fluctuations is equal to the area under the power spectrum, we get

∫_{allk_{z}} Ψ(k_{z})dk_{z} = ∑_{allk_{z}} < u^{2}(k_{z} >_{k_{z}band}≡< u^{2} >  (83)

Following Dewan, we make the assumption that the bandwidth is proportional to the central wavenumber

Δk_{z} ∝ k_{z}  (84)

and the assumption that all wavenumber components are limited by the previous saturation condition [eq 78] then we get

Ψ(k_{z})k_{z} ∝ u_{o}^{2}(k_{z})  (85)

and therefore

Ψ(k_{z}) ∝ \frac{4N^{2}}{k_{z}^{3}}  (86)

which tells us the predicted slope of the spectrum.

We can get to this result much more easily if we use a similitude argument. If we assume that the saturation due to the shear instability determines Ψ, then Ψ must depend on the critical shear S_{cr}^{2} = 4N^{2}.

The only other variable that Ψ could depend on is k_{z}. Using the fact that the dimensions of Ψ are [L^{3}T^{−2}] and those of S_{cr} and k_{z} are [T^{−1}] and [L^{−1}] then we immediately obtain

Ψ ∝ \frac{S_{cr}^{2}}{k_{z}^{3}} = \frac{4N^{2}}{k_{z}^{3}}  (87)

It is interesting that there seems to be less assumptions made in the latter derivation.

Below is a summary of the predictions made by Dewan for the frequency and wavelength spectral slopes.

Ψ_{u,v}(k_{z}) ∝ k_{z}^{−3} (m^{2}/s^{2})/(rad/m)  (88)
Ψ_{u,v}(k_{hor}) ∝ k_{hor}^{− \frac{5}{3}} (m^{2}/s^{2})/(rad/m)  (89)
Ψ_{u,v}(ω) ∝ ω^{−2} (m^{2}/s^{2})/(rad/s)  (90)

It is this latter equation that will be tested in the experimental results section below.

B. Doppler Spread Theory

The Doppler spread theory of gravity wave saturation depends naturally on the Doppler effect. In short, the Doppler effect describes the change in frequency of a wave when the observer is in motion relative to the medium. The frequency of the wave increases if the distance betwen the source and the observer decreases, while the frequency decreases if the distance between the observer and the source increases.

The Doppler shift of gravity waves is concerned with the phase velocity of the gravity wave and the mean flow. As mentioned earlier, a gravity wave will propogate upwards until it reaches the critical layer. At this point the phase velocity is equal to the mean flow, or from the point of view of Doppler spread theory, the frequency of the wave is Doppler shifted to zero.

This effect can easily be described by way of a diagram (Figure 16). The veritical axis is northward velocity and the ordinant is eastward velocity.

Fig. 16. Hines “critical circle”, [Hines, 1991]

As the gravity wave is propagated upward from the troposhere it grows in amplitude until it reaches a critical layer. This critical layer is represented by a circle on the diagram (Figure 16). The reason that it is a circle is because it is the component of the phase velicoty in the direction of the the mean flow that causes the Doppler shift. This means

V_{phase} = kV cos(θ).  (91)

For θ ≠ 0 the component of phase velocity in the direction of the mean flow is equal in manner to the mean flow, and the wave breaks.

The Doppler spread theory of Hines [1991] points out that the wave interactions with the mean flow alter the mean flow, and hence the critical circles are much broader than the thin lines shown in the figures (Figure 17). The shaded part is the perturbation due to the gravity waves themselves. Thus waves are eliminated at a greater rate than predicted by linear saturation theory.

Fig. 17. Hines “critical circle”

The Wienstock nonlinear cascade theory has been shown by Hines to be similar to Doppler Spread Theory, and therefore will not be examined here.

Turbulent Breakdown of Internal

Atmospheric Waves:

II Experimental Results

This section describes actual measurements made here at Delaware Observatory over the last few months. The hydroxyl airglow (an eight kilometer thick layer centered at a height of 87 km) is used as an indicator of gravity wave activity. The Power Spectrum of the fluctuations in the OH airglow is calculated and compared with the gravity wave theory discussed in the previous section.

VII. Power Spectrum Estimation

Of frequent concern to researchers is the amount of energy in a time series present at any particular range of spectral components. This falls under the heading of Power Spectrum Estimation. There is a veritable plethora of methods for estimating the power spectral density (PSD) of a time series. The most basic is the periodogram method. The power spectrum is defined as the squared modulus of the complex value of the spectrum calculated by the DFT. By the convolution theorem, this is equivalent to taking the autocorrelation of the time series and transforming it to the frequency domain.

There are difficulties in estimating the power spectral density (PSD) however. It is an unfortunate fact that there is always noise inherent in any real data gathering. To handle this, one calculates the significance of the spectral peaks with the assumption of a gaussian distribution of white noise (ie. a uniform spectrum). It turns out that the standard deviation of the time series is equal to the area under the power spectrum (Avery, 1993), where the standard deviation δ of a function h(t) is defined by the variance σ^{2}. In an attempt to reduce the variance, techniques such as spectral averaging are used. In this case, the time series is broken into K consecutive series of M points each such that K · M = N where N is the total number of points in the time series. A periodogram is calculated for each of the K series and these K PSD’s are summed together to give a final PSD of the time series. The variance is reduced by the square root of K (Avery, 1993), however this is at the cost of frequency resolution, since there are then only M frequencies sampled instead of the N frequencies.

VIII. The OH Layer

Chemical reactions in the atmosphere can store the radiative energy from the sun during the day, and this energy can be re-emitted at a later time. During the day, this airglow is largely swamped by the brightness of the scattered sunlight, but during the night, this chemiluminescence contributes significantly to the brightness of the night sky, especially in the near infra-red. The brightest source in the infrared is the vibrationally excited hydroxyl airglow (Chamberlain, 1961), called the Meinel Bands for the pioneering work of Meinel (1950) in observing the rotational-vibrational lines. The hydroxyl layer lies near the mesopause and has a thickness of about 8 km (Baker, 1988). This region of the atmosphere is too high to be measured in situ with balloons and such, and is too dense for a satellite’s orbit. Rockets can probe this region, but cannot gather extensive data. Remote sensing is the only way to take extensive data. In the main nighttime destruction mechanism for ozone, the layer is formed principally by the following reaction (Bates and Nicolet, 1950),

H + O_{3} → OH + O_{2} + 3.4eV  (92)

This exothermic reaction can excite vibrational level of hydroxyl up to the v′ = 9 level, corresponding to 3.3 eV. The v′ = 10 level ( 3.6 eV ) is not excited by this reaction and negligible emission is seen from this level. The atomic hydrogen is replenished in

OH + O → H + O_{2}  (93)

The (4-2) and (3-1) transitions are actually the best transitions to observe in the night airglow because are they easily detected. There have been numerous investigations of the hydroxyl airglow (for example Lowe 1969, Krassovsky and Sefov 1972, Baker and Stair 1988, Mcdade and Llewellyn 1988, Takahashi et al. 1992). The intensity of an individual line for an OH vibration-rotation band is given by (Herzberg, 1955)

I(v, v′, J, J′) = \frac{(2J + 1)N(v′)A_{v,v′,J,J′}e^{− \frac{hcF (J′,v′)}{kT}}}{Q(v, T)}  (94)

where J is the total angular momentum, N(v) is the population of the v vibrational level, A is the Einstein coefficient for the transition, h is Plank’s constant, c is the speed of light, F is the rotational term value, and Q is the rotational partition function. If the effect of a passing gravity wave on the OH layer is to raise or lower the layer adiabatically, then one can relate the height of the layer to temperature changes through the adiabatic lapse rate, which for a height of 85km is about 10^{o}/km. While the change in density is the most significant cause for a change in intensity, changes in the temperature will also affect the reaction rate coefficients in the processes that produce the layer, and therefore the intensity of radiated light will change. The increase in density also affects the reaction rates, and even the amount of hydroxyl within the field of view. Thus, measuring the changes in the intensity of the hydroxyl airglow gives a valuable insight to the effect of a passing gravity wave.

IX. Analysis of Delaware Observatory Data

The measured OH intensities used for this paper were taken with the University of Western Ontario’s SCanning Radiometer. It is a simple imager of the OH emissions. The instantaneous field of view is 1 km in diameter, and it takes a 16 by 16 point square of measurements in a raster scan every 60 seconds. The detector is a thermo-electrically cooled single element InGaAs (Indium Gallium Arsenide) device. It has a silicon window to filter high frequency light (cut-off at 1700 nm) and is limited at the low frequency end by the detector (1200 nm). The instruments is computer controlled, and runs automatically, turning itself on at a solar depression angle of 6^{o}. The intensity of the emissions have not been absolutely calibrated, thus the units are arbitrary, but proportional to Rayleighs.

Fig. 18. The intensity of the OH airglow as measured by UWOSCR for Sept 27, 1995

The signal to noise of the instrument was measured by covering the instrument and running it for a night. This also gave the DC level. The standard deviation was found to be ±10 arbitrary units [Huang, 1993]. The DC level also changed with the outside temperature, since the detector is cooled to a constant 40^{o} less than ambient. This could result in a systematic error in UWOSCR measurements, but since temperature is slowly varying, it can be removed from the signal in this analysis. All the time series used here were detrended with a least square fit to a second order polynomial.

They are taken only during the night hours, as the hydroxyl layer cannot be observed in daylight. There were also a few nights where poor weather prohibited observing. The time series were approximately 9 hours long, sampled every minute. An example is given in Figure 18.

The rising and setting of the sun contaminated the beginning and ending of the time series, and these were removed. The dates used in the calculation of the power spectra density range from September 1, 1995 to November 13, 1995.

Fig. 19. The average power spectrum for the nights listed above.

The PSD shown in figure 19 was the average of 29 nights. Each night was approximately 9 hours long, sampled every minute. The time series was detrended by removing 2nd-order polynomial fit, windowed with a hanning window and Fourier transformed to the spectral domain. Each night’s PSD was averaged to give the final result. The slope was calculated with a least squares fitting routine using the points ranging from 4 hours to 10 minutes in period. The slope found was -1.99 which is agreement with Dewan’s similitude prediction.

X. Conclusion

The recent data from Delaware Observatory confirms several key points in gravity wave theory. The spectral slope predicted by Dewan is indeed seen in the OH airglow fluctuations over the inertial range of 9 hours to 10 minutes. Also apparent on that graph is the limiting high frequency known as the Brunt Vaisala frequency. The power levels off to a noise floor at a frequency corresponding to a 5 minute period. This confirms the dispersion relation which implies that no gravity waves can exist with frequencies higher than maximum bouyancy frequency.

Clearly the turbulent dissipation of gravity waves is of fundamental importance in the understanding of earth’s atmosphere. Further work is required in understanding and exactly how the gravity wave field

affects the momentum and energy budgets of the atmosphere.

References

[1] Avery, S. (1993), The periodogram and its statistical properties, School on Time Series Analysis, 1993 C.E.D.A.R. Workshop.

[2] Baker, D.J., and A.T. Stair (1988), Rocket measurements of the altitude distributions of the hydroxyl airglow , Physical Scripta, Vol.37, Gn-622

[3] Beer, T. (1974), Atmospheric waves, John Wiley and Sons, New York

[4] Chamberlain, J. W. (1961), Physics of the Aurora and Airglow, Academic Press, New York, U.S.A.

[5] Dewan, E.M., and R.E. Good (1986), “ Saturation and the ‘Universal’ Spectrum for vertical profiles of horizontal scalar winds in the atmosphere”, J. Geo. Res., vol. 91, no.D2, 2742-2748

[6] Fritts, D.C., and P.K. Rastogi (1984), “Convective and dynamical instabilities due to gravity wave motions in the lower and middle atmosphere: Theory and observations” Radio Science, Vol.20, no. 6, 1247-1277

[7] Hines, C.O. (1960), “Internal atmospheric gravity waves at ionospheric heights”, Canadian Journal of Physics, vol.38, 1441-1481

[8] Hines, C.O. (1967), “On the propagation of atmospheric gravity waves through regions of wind shear”, J. Geophys. Res., vol.72, 1015-1034

[9] Hines, C.O. (1974), “The upper atmosphere in motion”, American Geophysical Union

[10] Hines, C.O. (1991), “The saturation of gravity waves in the middle atmosphere. Part II Developement of Doppler-Spread Theory”, J. Atmo. Sci.,vol 48, no. 11

[11] Garret, C. and W. Munk (1975), “Space-time scale studies of internal waves: A progress report”, J. Geophys. Res., vol.80, no.3, 291-297.

[12] Huang, Y. (1993), Analysis of gravity wave activity in hydroxyl airglow, M. Sc. Thesis, University of Western Ontario, London, Canada

[13] Lewellen, W.S. (1970), “Use of invariate modelling”, Handbook of Turbulence, Walters Frost and Trevor H Moulden, Eds vol.1, Plenum, 237-280

[14] Meinel, A.B. (1950), Astrophys. J. 111, 434-444

[15] Patnaik, P.C., F.S. Sherman and G.M. Corcos (1976), “A numerical simulation of Kelvin-Helmholtz waves of finite amplitude”,J. Fluid Mech., Vol.73, 215-240

[16] Peltier, W.R., J.R. Halle and T.L. Clark (1978), “Evolution of finite amplitude Kelvin-Helmholtz billows”, Geophys. Astrophys. Fluid Dyn., vol.10, 53-87

[17] Phillips, O.M. (1977), The dynamics of the upper ocean, Cambridge University Press, New York

[18] Swarztraube ,PN., and R.A. Sweet (1979), “Efficient FORTRAN subprograms for the solution of seperable elliptic partial differential equations”, Assoc. Comp. Mach., Trans. Math Software, col.5, 352-364

[19] Sica R. (1995), Meteorology Lecture Notes, The University of Western Ontario

[20] Sykes, R. I. (1982), “A numerical study of breaking Kelvin-Helmholtz billows using a Reynolds-stress turbulence closure model”, J. Atmo. Sci., vol 39, 1506-1520

[21] Thorpe, S.A. (1973), “Experiments on instability and turbulence in a stratified shear flow”, J. Fluid Mech., Vol.61, 731-751