A.R.A.P. REPORT NO. 169
AN INITIAL TEST OP THE
APPLICABILITY OP INVARIANT MODELING
METHODS TO ATMOSPHERIC BOUNDARY
LAYER DIFFUSION
Coleman duP. Donaldson
and
Glenn R. Hilst
October 1971
to
Environmental Protection Agency
National Environmental Research Center
Division of Meteorology
Research Triangle Park, North Carolina 27711
EPA Contract No. 68-02-0014
AERONAUTICAL RESEARCH ASSOCIATES OF PRINCETON, INC
50 Washington Road, Princeton, New Jersey 085^0
-------
TABLE OP CONTENTS
Page No
Nomenclature 111
1. Introduction 1
2. The Basic Diffusion Equations and Method of
Closure 4
2.1 Closure at Higher-Order Moments 7
2.2 The Method of Closure at Second-Order
Moments 8
3. The Invariant Diffusion Model 12
4. Some Relationships Between Second-Order
and First-Order Closure Models 18
4.1 K* from First- and Second-Order Closure
mModels in a Neutrally Stratified, Constant
Shear Stress Layer 22
5. Initial Tests of the Model 26
5.1 An Example' 27
5.2 An Initial Test of the Effect of Stability 29
5.3 The -Sensitivity of Diffusion to the Scale
Length 31
6. A Calculation of the Time Required for Turbulence
to Become Fully Developed 35
7. Conclusions • 37
8. References 39
Figures 1 through 12 40
Appendix A 52
ii
-------
Nomenclature
a,b constants
C, C pollutant concentration (mass fraction)
c specific heat at constant pressure
g acceleration of gravity
i subscript index
222
K u'+v'+w'
K*, K*, K*, K* eddy diffusivities for heat, momentum, and
H m y z matter
k von Karman's constant
L Monin-Obukhov length
n an integer, 1, 2, 3,...
p pressure
Re. turbulent Reynolds number
R. Richardson number
S source or sink strength of pollutant (per unit
volume and time)
T, 1" mean and fluctuating values of potential
temperature
t time
u, v, w time-averaged components of fluid motion
u1, v1, w1 turbulent components of fluid motion
x, y, z Cartesian coordinate system
JSr molecular diffusivity
a, P Monin-Obukhov coefficients
A, A. diffusivity and dissipation scale lengths
p. viscosity
p,p second viscosity coefficient
p, p fluid density
An overbar indicates a time-averaged value of the quantity; a
prime indicates an instantaneous fluctuation of that quantity
about its average value. Averages of products of fluctuations
of two or more quantities are sometimes referred to as
correlation terms.
d/dt local time derivative
d/dx-j_ local space derivative
d d 3 S S
-TT- = T—- + U -=r— + V T-— + W ^r—
at at dx dy dz
iii
-------
1. INTRODUCTION.
Recent emphases on the alleviation and prevention of
significant'air pollution problems have placed a new premium
on reliable methods, for predicting the transport, diffusion,
chemical reaction, and deposition of airborne materials
emanating from man-made sources. Since the capacity of the
lower atmosphere to. disperse airborne wastes to harmless
concentrations varies through two or more orders of magnitude
within a single day and from one place to another, this
variable alone determines to a large extent the feasible
utilization of the atmosphere for waste disposal.
A variety of techniques for predicting atmospheric
dispersion have been brought forward to address this waste
management problem. These techniques vary in sophistication
from simple kinematic descriptions of the accumulation .of
pollutants.within a hypothetically closed box to highly
complicated computer solutions of -the combined equations of
atmospheric motion and complex photochemical chain reactions.
At the present stage of development, all of these models
require a large content of empirical knowledge, a feature
which makes their application to new meteorological conditions,
different'terrain, or .new pollutant .source geometries highly
questionable.
The work reported here represents an initial attempt to
remove some of this empiricism from atmospheric diffusion.
calculations by an appeal to a more rigorous physical descrip-
tion of the process of turbulent diffusion in the lower
atmosphere. In particular, it has appeared worthwhile to
investigate the applicability of the method of invariant
modeling of turbulent fluid flows to the lower atmosphere and
its diffusion capabilities. These techniques, which have been
successfully applied to classical shear flows (such as the free
jet, the plane shear layer, and the free vortex) depend upon
-------
a closure of the governing equations at the level of the
equations for the second-order correlations or moments of
turbulent fluctuations. These moments contain, among others,
the correlations responsible for the transport of momentum,
heat, and-matter by the turbulent motion. The success of these
techniques in predicting the structure .of turbulence in classi-
cal shear flows with rather a wide variety of geometries, e.g.,
from free jet to free vortex, suggested that the modeling might
be profitably carried over to the more complicated atmospheric
flows and the turbulent diffusion of a pollutant.
Discussions of this possibility with the Division of Meteor-
ology of the Environmental Protection Agency led to the conclu-
sion that a1test of the invariant modeling technique, in its
simplest form, should'precede any further significant develor
ment of these concepts for atmospheric problems. It was,
therefore, agreed that a modest effort to test the model for
the case of an infinite cross-wind line source in the lower
atmosphere would be undertaken as a first demonstration of
the general utility of this technique. Successful demonstration
of the model for relatively simple source and atmospheric flow
conditions would presumably justify additional refinement and
development of the technique. In particular, the ability of
the method to provide more precise physical insights into the
complicated mechanisms of turbulent diffusion in .arbitrarily
stratified atmospheres was to be pursued at a later stage.
The initial tests of the earliest model have been
completed and are reported here. These have been successful
beyond reasonable expectations. In particular, using only the
mean wind and temperature profiles for the lower atmosphere and
a predetermined scale length, the vertical dispersion of a
pollutant has been predicted well within a factor of two of
observed values out to distances of about two miles. The
predictions were somewhat more exact for neutral and unstable
temperature stratification than for stably stratified atmos-
pheres, but even the latter were remarkable considering the
-------
general lack of any satisfactory theory for diffusion in
stably stratified atmospheres.
Even with these successes, the present effort has been
most instructive. A greater understanding of the physics of
turbulent transport has been achieved and has led directly
to further improvement and refinement of the model. The
dependence of the scale on stability.has been clarified and,
although this dependence is not fully.understood, it is clear
that the ratio of the scale of the breadth of a pollutant
plume to the integral scale of the velocity fluctuations in
the atmosphere plays a very important role in the physics-of
the dispersal process.
With the successful demonstration of this primitive
form of the invariant model, attention has been turned to the
refinement of the model in the light of recent numerical
studies of classical shear flow, the construction of a three-
dimensional diffusion model, and the application of these
results to initial calculations of diffusion in arbitrarily.
stratified atmospheres. Among the refinements planned is the
inclusion of the coriolis forces in the equations of motion so
that diffusive effects within the Ekman layer may be calculated.
-------
2. THE BASIC DIFFUSION EQUATIONS AND METHOD OF CLOSURE
Since the advent of the .Navier-Stokes equations and
Reynolds' introduction into these equations of the mean and
fluctuating components of the vector and scalar quantities
involved, various attempts have been made to construct
predictive models for the structure and diffusive character-
istics of turbulent fluid flows. The unclosed character of
these equations and the lack of powerful computers has, until
recently, dictated the,use of only the very simplest schemes
for>a modeling closure of the equations of motion. These
methods were not without merit and fruitful results; to the
contrary, much of our present engineering competence in areas
such as fluid.flow in pipes, boundary layer flow, and even
atmospheric diffusion processes has depended heavily upon
these early insights and the organized empiricism to which
they gave rise.
Although it is not our intention to review the develop-
ment of methods for computing turbulent shear flow (the
reader-interested in such a review should consult Reference 1),
it is instructive to recapitulate the most frequently used
method for closure of the governing equations as background
to the present method of closure developed by Donaldson and
termed, "invariant modeling." The earlier methods of closure
are readily illustrated by manipulation of the differential
equation for the passive dispersion of matter in an incom-
pressible fluid undergoing turbulent motio.n.
Let C be the mass fraction of a suspended gas or
aerosol, p the.density of the fluid, and let us adopt a
cartesian coordinate system (x,y,z) with the components of
fluid flow u, v, w along each of these directions. Then the
governing equation'for C , following the motion of a tagged
fluid element, is
-------
dC SC 0-C BC o-C
eft = P Bl + Pu to + PV * + Pw 5¥
where
-------
The basic diffusion equation which it is customary to
use is obtained from Eq. (4) by noting the experimental fact
that gradients of mean quantities in directions normal to
the mean motion are generally much larger than gradients in
the direction of this motion, so that the term
in (3) may be neglected. Under this assumption, which we
shall adopt only for convenience since it is not required, the
diffusion equation that one wishes to solve is
This equation may be solved, subject to appropriate initial
and boundary conditions, if we specify some relation between
the second-order correlations C'v' and C'w1 and the mean
distributions of concentration and velocity. Such a simple
closure is useful, provided the flow is such that the. dynamics
of the second-order correlations can readily follow or "track"
changes in the mean quantities.
Arguing from the physical rationale of a "mixing length"
analogous to a molecular mean free path, Prandtl provided an
early closure scheme by writing
and
, -PK* (6)
= -pK* (7)
where K* and K* are eddy dif fusivities , quite analogous
to but generally much larger than the molecular diffusivity
& . Substitution of Eqs . (6) and (7) into (5) yields the
well-known diffusion' equation
-------
Strictly speaking, Eqs . (6) and (7) only define K*
«y
and K* and closure is not achieved until these parameters
Z
are specified as a function of position and the fluid flow
and concentration fields. A very large literature devoted
to these specifications of eddy diffusivity continues to
this day.
At the present, our concern is not with the successes
or failures of this method of closure but, rather, with clear
identification of the method itself. Stated generally, the
principle involved is that closure may be.achieved by modeling
the nth order correlations in terms .of the lower order
correlations. Prandtl's method, which has just been described,
invokes this principle at its lowest possible level, i.e.,
n = 2 . Thus, the prevalent methods of closure retain only the
minimum information inherent in the governing equations, since
any effects of second- or higher-order moments are assumed to
be explicitly contained in the mean values and the scale lengths
of the flow and concentration fields. Since n = 2 is the
minimum level at which closure may be obtained, we shall term
methods such as Prandtl's "first-order closure."
2.1 Closure at Higher-Order Moments
Since, in principle, closure may be achieved at any order
of correlation, the possibility of retaining significant
information by closure at the third- or higher-order moments
appears worth investigation. The techniques for deriving
equations for the higher-order moments have been known for
some time (for a complete derivation, see Reference 1, p. 250,
et seq) . To illustrate the technique, we note, for example,
-------
= c
and
dv' _ iy_ _ li fir}}
at ~ at at (W)
__
at
From the equation-of motion for v and the diffusion
equation for C , we form the mean values 5v/dt and
subtract these from the original equations to find av
and ac'/at, and then multiply by C' and v', respectively,
and add. We then take the time average of this equation to
find aC'v'/at . When this process is carried out, we always
find the resultant equation contains the next higher correla-
tion or moment, e.g., v'v'C' ; hence, with each prediction
equation for the nth correlation, we invoke a new unknown,
the n+1 correlation term. Closure, therefore, depends upon
successful modeling of the n+1 correlation terms in terms of
the lower correlations for mean values of the flow field. It
is to this problem that the methods of invariant modeling have
been primarily addressed.
2.2 The Method of Closure at Second-Order Moments
Against the background considerations outlined above, we
may note generally that any attempt to model turbulent fluid
motions and the diffusion of a passive pollutant involves
consideration of seven variables:
u, v, w the components of fluid motion
T the fluid'temperature
p the pressure
p the fluid density
C the concentration of pollutant
-------
To construct a model we must have seven governing equations
for these variables. These are supplied by the Navier-Stokes
equations (three equations), the diffusion equations for heat
and matter (two equations), the equation of mass continuity,
and the equation of state. The latter equation relates p ,
T , and p and reduces the independent variables to six.
Prom.these governing equations, we may develop the equations
for any higher-order moments that we wish to consider. If we
wish to attempt closure at the second-order moments, i.e., model
any third-order moments in terms of second- and first-order
moments, and'we have six independent variables, we must, in
general, develop equations for the following twenty-seven first-
and second-order moments.
First-order moments:
u, v, w, p, T, C 6 equations
S'econd-order moments :
u' ,V , w' , T' , p' , C'* 6
u'v', v'w', u'w' 3
v'T', w'T' 3
u'p', v'p'. w'p'. p'T'
u'C', v'C', w'C', C'T', C'p' _J5
Total Equations 27
In principle, and given successful modeling of third-order
moments in terms of lower order moments, we should be prepared
to contend with 27 simultaneous partial differential equations,
a prodigious task indeed.
Fortunately, the required equations are reduced signifi-
cantly by several considerations. First, two of these moments
2
do not appear in the derived moment equations, namely, p' and
2
C1 . Equations may be'written for these variances, but they
do not feed back and are, therefore, not required for closure.
Second, we may eliminate the following moments by assuming a
-------
10
horizontal straight-line (nondivergent) flow and an
incompressible fluid:
v, w, u'v1 , v'w' , v'T'
a further reduction by five equations.
As a further step in reducing the number of equations,
p has been eliminated by linearizing the equations, using
essentially the Boussinesq approximation for small depart-
ures of p , T , and p from their values in an adiabatic
atmosphere at rest. According to the methods of invariant
modeling, the covariances of pressure fluctuation and each of
the other five variables are modeled in terms of the remaining
moments. This procedure which is described in the appendix
removes u'p', v'p', w'p', p'T', and p'C' from separate
consideration.
These considerations reduce the number of simultaneous
partial differential equations from 27 to 14. In the early
version of the invariant model, one more equation has been
eliminated by retention of the assumption of a continuous
pollutant source and that du'C'/dx is negligible compared
with dv' C' / c)y and Sw'C'/dz in this case. Under this
assumption, the turbulence model developed from second-order
closure involves 13 simultaneous partial differential equations
for first- and second-order moments. In addition, this method
of closure invokes two scale lengths: the scale of the turbu-
lence A which we treat as a free parameter, and the dissipa-
tion scale length X which is related to A by an algebraic
2 2
expression involving two constants a and b , u' + v' +
2
w' = K , the fluid density p , and the molecular viscosity
I_LQ (see page 17 ) .
In the general case of an active pollutant, e.g., a
pollutant which is hot with respect to the fluid temperature,
the 13 partial differential equations must be solved
simultaneously. However, for passive pollutants there is no
feedback from the pollutant diffusion equations to the fluid
-------
11
flow equations. Under this assumption, the equations for C",
v'C', w'C', and C'T1 may be decoupled and solved separately
after solutions for the remaining nine moments associated with
the fluid flow and temperature have been obtained. In the
early tests of the model reported here, a further economy
of computer operation has been realized by taking known
values of u(z) and T(z) as fixed inputs to the program.
Under this mode of operation, the equations for u and T
are not solved, reducing the number of simultaneous partial
differential equations for the turbulence component of the
model from 9 to 7.
Beyond the general considerations discussed above, the
development of the invariant model.for turbulence and
diffusion has depended upon techniques for modeling third-
order moments in terms of first- and second-order moments.
The development of these techniques, has been discussed by
Donaldson and his colleagues [2, 3, 4, and 5] (Ref. 5 is
appended to this report) and will not be repeated in detail
here. (A more rigorous report of the development of these
techniques and their application is now planned, following
this successful initial test of the concepts involved.)
-------
12
3. THE INVARIANT DIFFUSION MODEL
When the procedures of invariant modeling are system-
atically applied to the equations of motion for the atmospheric
boundary layer and the diffusion equations for a (scalar)
contaminant, the simplest set of simultaneous partial differ-
ential equations for the prediction of the diffusion of a
contaminant emitted continuously from a point source is as
follows:
1. The diffusion equation*
dC_ \ . / dc" \
a_ U _P_ _ p c'w'
5 \ I
(12)
u -r-K- = S + ^- p. ^ - p C ' v ' + ~ u ^ - p C ' w '
>o d ro p I £z I o 3Z Fo p
This equation states that the rate of change of the average
concentration C , following the motion, is equal to any local
source or sink plus the local diffusion due to molecular motions
2
^ _ "C /dy and the local flux divergence -S(p C'v'
^o p' ro p '
(We retain the molecular terms throughout this development so
that the equations may be solved right down to the solid
boundary .) \
2. The turbulent flux equations
PV
3 ^7 (PO
Av/K _±L_ ) + |- ' - *-*—*-
' oz / oy
BC'w1 \ p^^K
» f p AV/K ^J2— - ^
3y2 Sz2
2C'v
JB— ] (13)
*This development is made under the assumption that the .
Schmidt number p £r/\L = 1, or u. = p £r, and that p = p (z)
only. oo oo o o
-------
13
and
N / dc' v'
4- [ p Av/K T^—
o dz
X /
i 1 o A /f^
/ a2cTwT
1 LJ. 1 P
"° I ay2
- K + ^2)
pns
H- ^- C ' T '
Lo p
BCDV' \ PO
k/ 1 W
^^Q , w ,
P
' 3,2 "
v^C
f ' W '
A CpW
2C'w' \
P \
x2
" 1 3\ 2 /3Po\2lc^
Lp^^--^(5rj j p
(14)
These two equations begin to provide some insight into
the physics of turbulent flux of a contaminant. The first
term on the right,of Eq. (14), -pow'w'dC /5z , is a produc-
tion term for C'w' which shows this flux is generated out
P 2
of the interaction of the turbulent motion w1 and the mean
gradient of concentration. The only other possible production
is the last term of (14), pogC'T'/To , which is a production
term if C ?T' is positively correlated, i.e., positive fluc-
tuations of C are associated with positive temperature
fluctuations, and is a dissipation term if C'T1 is negative.
This term, of course, simply reflects the effects of buoyancy.
The 'second and third terms on the right side of (13) and (14)
represent the diffusion of the flux along the direction of the
flux. The fourth and fifth terms represent the transformation
of flux in one direction into flux in the orthogonal direction.
The sixth term represents the loss of flux due to the tendency-
to-isotropy of the turbulent flow. The seventh term represents
the molecular diffusion of the flux and the reduction of the
-------
14
magnitude of the flux by the effects of molecular motion.
The eighth term in (14) reflects the fact that, in the presence
of a density gradient, the divergence of the flow is nonzero
and leads to a small production of flux in the vertical direction,
These equations have produced one new term, C'T1 ; hence
the following equation for C'T' is required.
3. The pollutant-temperature correlation.
i I m i ?sfi / ^p i m
iL . -p ^F S> - p 5^ |T 3 f AyK ±£^
c ro dz ro p ^z ^y \ ^o ^y
2C'T'
"^~
(15)
Again, the first two terms are production terms for C'T' ,
the second two are eddy diffusion terms, and the fifth repre-
sents molecular diffusion and a reduction of C'T' .
These four equations provide the diffusion model for
material released continuously from a point source. Since
the present test of the model is for an infinite cross-wind
line source, we may simply note that all of the terms in
Eq. (13) are zero for this source configuration. Equations
(12), (14), and (15) then constitute the diffusion model.*
In order to solve these three equations, we require prediction
equations for u(z), T(z), u'2, v'2, w'2, u'w', u'T', w'T' ,
2
and T' . Since these are the same equations required for
prediction of the turbulent flow, independent of the presence
of a passive pollutant, we may utilize the invariant model of
*We have named this diffusion model LPD, the acronym for Line
Pollution Diffusion. The model for a point source has been
designated as SPD for Spatial Pollution Diffusion
-------
15
turbulent boundary layer flow directly. This model has been
discussed extensively by Donaldson, including the reprint
attached as Appendix A to this report. (For another discussion
of the physical meaning of the terms appearing in this model,
the reader is referred to Lumley and Panofsky [6].) For
completeness, we repeat the turbulent model equations as used
for the present tests.
4. The equations for the mean wind and temperature.
PC H -^o 0 - Iz ^7 <16>
5. The turbulent energy equations
o
K
3
-------
dt
apr
' w
oz
+
o
w'w'
O
2
PO
6. The momentum and heat flux equations
du
'
A
+
"Sz
u'w1 o
o
o
u'w1 . °
PO \
/dp
Si
16
T
w 'T
o
(20)
p K
§-u'T'
(21)
o
'1
p A/K u'T
+
o
(22)
-------
dw ' T" — ; — : dT _ d / . « dw ' T '
= - w'w- -- +
O
oz I ° dz
A
7. The variance of the temperature fluctuations
PnTT— = -2p0w'T' T + /
T'2
O
17
Following Rotta [7], we complete the model by writing
A = X/a + bReA (25)
where a and b are constants, and
p
A
For our present purposes, A has been retained as a free
parameter.
(26)
-------
18
SOME RELATIONSHIPS BETWEEN SECOND-ORDER
AND FIRST-ORDER CLOSURE MODELS
One of the purposes of second-order closure Is to
retain more specific information regarding the structure
of turbulent flows, out of which, hopefully, increased
physical insights can be gained for the improvement of
first-order closure models of .turbulent fluid flow. These
relationships between first- and second-order closure models
are most readily,seen by comparing the predictions of the
two types of models for specific physical processes. For
our present purposes, we choose to compare predictions of
the vertical flux of-.horizontal momentum p u.'.w' and' the
vertical flux of heat • p c w'T' which are given by
p u'w' = - p K
o o m
and
(27)
(28)
where T is potential temperature and c is the specific
heat at constant'pressure. Equations (27) and (28) are best
interpreted as the definitions of K* , the eddy diffusivity
for momentum, and K* , the eddy diffusivity for heat, both
of which may be functions of height, stability, wind speed,
surface ro'ughness, and each other. We may, therefore, focus
attention on comparisons of first- and second-order closure
models for K* and K* ; the latter are defined as
171 rt
K* = - u'™' (29)
m (au/dzl
and
K| . - _pl_ (30)
H
-------
Since all of the terms on the right-hand side of Eqs.
(29) and (30) can be, and have been, measured in the lower
atmosphere, we have relatively good empirical evidence of
the values of K* and K* in that layer and their general
m n
behavior as a function of height, wind speed, stability,
and surface roughness, The predictions of first-order and
second-order closure models for K* and K* may, therefore,
be tested against reality, as well as against each other.
Various empirical and semi-empirical forms for K*' and
K* have been proposed (cf., Priestley [8], Pasquill [9]>
n.
Lumley and Panofsky [6], and Pandolfo [10]). Prom these, we
choose to compare the following forms of first-order closure
models:
"Profile" Model
Km = k2*2 H (1? aRi)±2 (31)
111 \J £i \ -L./
and
"Similarity Theory" Model
K* - kcz* & (I + £2,rc (32)
where R. is the Richardson number and L is the Monin-
Obukhov length. These parameters are defined by
and
, (33)
1 T
L = -
kg w'T1
The choice of sign in Eq . (3D is determined by the sign of
the Richardson number ; the upper sign is chosen when R- >_ 0,
and the lower sign when -0.05< R. < 0 . (Since both Eqs.
(3D and (32) are restricted to a relatively small range of
-------
20
stability about'neutral (R. = 0, L = oo), we shall restrict
the comparison between first- and second-order closure models
to that range, also.)
The predictions of K* and K* from the invariant
m H
model may be derived directly from Eqs. (21) and (23), if we
assume an equilibrium condition has been reached, i.e.,
Su'w'/3t and ow'T'/St are equal to zero. Neglecting all
of the molecular terms except the molecular dissipation
(which is tantamount to restricting the comparison to heights
well above the laminar sublayer), we solve for u'w1 and
vTnFr from Eqs. (21) and (23) and form K£ and K* . The
results are:
"Invariant" Model
K» = 2 ! 1 2 (35)
'*o_ ^\dg
Po^+ A ; ^
and
p AT"1 Q ^ / __ /^-ui I T1 I \ rp P
w . ^ ££ - ^L 2—1 pA,/K ™ + —|r- T ' *
O7 P 017 I ^ OZ / P
C| = 2__\ . / o^ (36)
Focusing attention on Eq. (35) which we wish to compare
with Eqs. (31) and (32), we may note immediately that the
eddy diffusivity derived from second-order closure will
depend upon the following processes and coefficients:
1. The local^production of shear stress by the term
o _
w'
-------
21
The !Local diffusion of shear stress, by the term
The local production (or dissipation) of shear
stress, by the longitudinal turbulent heat flux as
measured by gu'T'/T (this is a production term
if u'T' is negative and is dissipative if u'T'
is positive) .
The coefficient of the essentially local dissipation
p
of shear stress by molecular effects, 2u. /p X • *
The coefficient of the tendency to isotropy of the
shear stress \/K/A . This coefficient depends upon
the square root of twice the local turbulent kinetic
energy which is measured by K , and. by the spatial
scale of the motion A which is not determined
I oc ally.
6. The local mean wind speed shear
Similar functional terms appear in Eq . (36). One important
distinction arises in the prediction of K* , however, in that
2
the buoyancy term (g/T )T" is always positive definite (or
zero) . The presence of temperature fluctuations always
enhances the value of K£ while they may increase or decrease
K* depending upon whether u'T1 is negative or positive. We
should also note that K* is not defined for du/Bz = 0 and
K* is not defined for T/^Z = 0. The latter condition is not
n
significant, since all of the terms in Eq . (36) are zero for
Sf/Bz = 0; i.e., there is no turbulent vertical heat flux in a
neutrally stratified atmosphere. This is not the case for K*
in i
when Su/dz =0, however, since only the local production of
shear stress necessarily goes to zero under this condition.
Shear stress produced elsewhere may be diffused into this region^
*X is of the order of 0.1 meters in the lower atmosphere but
may be as large as 10 meters in the stratosphere.
-------
22
and shear stress may be dissipated or produced in the
presence of fluctuations of the potential temperature and
the longitudinal component of velocity.
These portrayals of the complicated physical processes
which determine the effective eddy diffusivity point imme-
diately to the tremendous requirements which first-order
closure imposes upon either the Richardson number or the
Monin-Obukhov scale length. When viewed from the perspective
of second-order-closure, the surprising feature is that first-
order closure models have'been as successful as they have been!
Since the first-order closure models do not separate the
physical processes of production,, diffusion, dissipation, and
the tendency-to-isotropy of shear stress or heat flux, term-by-
term comparison is not possible, and we must depend upon
calculations and observations of K^ and K* to proceed with
any general comparison. However, at the present stage of
development, we may investigate the very special case of a
constant shear stress'layer in a neutrally stratified atmos-
phere, a situation for which K* is not defined. This is the
case for which the first-order closure models are most reliable,
and equivalence between first- and second-order closure
models can be expected.
4.1 K* from First- and Second-Order Closure Models in a
m
Neutrally Stratified, Constant Shear Stress Layer
Under the conditions ST/5z = 0 and du'w'/dz = 0
the form of K* from the invariant model is
m
and both the "profile" model [Eq. (3D] and the "similarity
theory" model [Eq. (32)] reduce to the form
-------
23
* = (38)
Equation (37) may be simplified further by noting that for
2
heights well above the laminar sublayer 2u. /p X «/K/A .
Then the invariant model predicts K* as
m
(39)
and we have equivalence among the "profile" model, the
"similarity theory" model, and the "invariant" model, if
2 2
k z
dz
w'2A
v/K
(40)
A physical interpretation of Eq. (40) is seriously
hampered by the obscuration of the physical construct of the
left-hand side and by the lack of specificity for A in
terms of fluid flow parameters, on-the right-hand side. A
p 2 _ p
simple dimensional comparison suggests that w' ~z (Su/oz)
-i *") O
and A//K ~(du/Bz)~ , but how k z du/^z is apportioned
between these two terms is not immediately evident.
From another point of view, the requirement that the
shear stress be constant with height leads to
o
= constant (4l)
2
If w' //K is constant, i.e., the vertical component of the
turbulent kinetic energy is a fixed fraction of the total
turbulent kinetic energy, we recover the logarithmic wind
profile (appropriate to this case) only if A~ z .
Since all of the terms in Eq. (40), other than A ,
are known or measurable, we can determine the value of A/z
t~\
from observations, of 3u/^z , W , and K in neutrally
stratified, constant shear stress layers of the atmosphere.
-------
Limited observations to date suggest A - 0.7z, but more
general observations over a variety of terrain are required
before this relationship can be confidently assessed from
experimental data.
At the present stage of development, it would, therefore,
appear most profitable to view Eq. (40) as the definition
for a reference value of A , say A , where A is inter-
preted as the value of this scale length which assures
equivalence between first- and second-order closure models
for the special case of a neutrally stratified atmospheric
layer for which the shear stress is constant.
We must note at this point, however, that in this early
version of the invariant model we have assumed there is a
single scale length A which is appropriate to both the
diffusion of shear stress and to the tendency to isotropy of
the shear stress. This assumption is not necessarily valid
and we should point out that the value of A appropriate to
Eq. (41) is the value of the scale length which enters into
the determination of tendency to isotropy. Modeling of
laboratory scale flows seems to suggest that this scale
length is about ten times the scale length appropriate to the
diffusion of shear stress. Further study of these two scale
lengths for atmospheric flow is required.
Beyond,this limiting comparison between first- and
second-order closure models, present attention is. being
directed to calculations, with the invariant model, which
match observed values of u'w' and w'T' as a function of
height in diabatic boundary layers. From such calculations
all of the inputs necessary to establish the vertical profiles
of the Monin-Obukhov length L and of a and p and their
dependence upon stability, wind velocity and shear, etc., can
be established. Once the credibility of the second-order
closure model has been established, it'will be an extremely
valuable tool in assessing the adequacy of first-order closure
-------
models. If these should prove adequate, particularly in
arbitrarily stratified atmospheres, the invariant model will
once again be useful in specifying the coefficients of first-
order closure models for these situations.
-------
26
5. INITIAL TESTS OP THE MODEL
Since the basic information normally available for
the atmospheric surface layer is the mean profiles of wind
and temperature, the method of solution adopted for the
calculations of turbulence structure was to decouple Eqs .
(16) and (17) and solve Eqs. (18) through (25) subject to the
observed (constant) values of u(z) and T(z). A small
222
initial turbulence (u' , v' , w' ) was introduced, and
this portion of the model was run until all the turbulent
energies and turbulent fluxes came into equilibrium; i.e.,
the local time derivatives went to zero.
In an earlier phase of this work,* data were obtained
from the Air Force Cambridge Research Laboratories (AFCRL)
which included mean wind and temperature profiles, as well
as careful measurements of the turbulent energies and fluxes
in the first 100 feet above a flat surface in Kansas. Although
no diffusion data are available in connection with the AFCRL
measurements, one of these runs was selected for diffusion
calculations. This run has been designated "Cote 40." It
was characterized by an unstable lapse rate from the surface'
to about 30 feet and was neutral above that height. The input
information u, T, and A and the calculated values of w'T'
and w'w' for this test case are shown in Figure 1.
In order to test the model under other stability condi-
tions, a second case was chosen for which the lower atmosphere
was stably stratified and for which detailed measurements of
the vertical distribution of a tracer material were available
out to 5000 feet from an elevated point source. These measure-
ments were made by Hilst and Simpson at Hanford and were
reported in the Journal of Meteorology [11]. Only the mean
wind and temperature profiles were available for this case.
*A.R.A.P. Report No. 162, NASA CR-111962, NASA Contract No.
NAS1-10192
-------
27
The input information and the calculated turbulence and
flux structure are shown in Figure 2.
With these choices, the initial tests of the model
bracket the stable to neutral or slightly unstable range.
Extension into the unstable case does not appear warranted
until a better match between observed and calculated turbu-
lence structures is achieved by the refinements presently
underway.
Having modeled the structure of turbulence, the diffu-
sion portion of the model was run and the predicted
concentration profiles and the moments of these profiles
were compared with experimental data. For this portion of
the test, the line source was placed at 200 feet above ground
level and the initial vertical concentration profile was
specified as a full gaussian with unit variance and a
central value of 1.0. Printouts of the predicted concentra-
tion profiles, the moments of these profiles, and the vertical
fluxes of material were required at distances of 500, 1000,
2500, 5000, and 10,000 feet from the source.
With this wealth of detail, the calculated results far
outstrip the available information from observational data.
However, a comparison with several experimental measurements
can be made by way of the variance of the concentration
2
distribution a (or a ) which is frequently reported. We
z z
can, therefore, readily check the model's order-of-magnitude
estimates of vertical diffusion. Beyond this, we must depend
upon a few observations and intuition in checking the predicted
concentration profiles directly.
5.1 An Example
In order to get an initial estimate of the behavior of
the diffusion model and to demonstrate the output information,
the first calculation was run under the assumption that a
-------
28
single scale length A was appropriate to the diffusion
simulation as well as for the generation of the field of
turbulence. Using the unstable to neutral case (Cote 40,
A =60 feet), the vertical profiles of concentration were
calculated to a distance of 10,000 feet downwind from the
source line. The basic results of this trial calculation are
shown in Figure 3 in the form of the concentration profiles
at x <= 500, 5000, and 10,000 feet; the values of the axial
concentration, Cn ; and the first four moments of the
^max
concentration distributions.
Inspection of the concentration profiles shows the
expected vertical transfer upward and downward, reflection at
the ground surface, and the establishment of an essentially
uniformly mixed layer below the height of the source. An
immediately disturbing feature, however, is the presence of
a narrow band of high concentration at the source height, a
feature which maintains abnormally high values of the axial
concentration. It would appear that with this initial choice
of A and the- constants a and b which relate A
lUcLX ITlciX
to the dissipation'scale length X , the model exhibits a
pathological tendency to minimize the turbulent flux near the
axis of the cloud, even though relatively large fluxes are
calculated outside this region. This problem, which is
discussed further in Section 5-3} has been of primary concern
and has led to a generalization of the scale lengths along two
lines: 1) the adoption of separate scale lengths for turbu-
lence generation and for diffusion, and 2) a reexamination of
the relative magnitudes of the scale lengths associated with
the dissipation of turbulent energy and fluxes as compared
with the pressure diffusion and tendency-towards-isotropy
terms. The first approach, which is largely empirical and is
aimed at an estimation of the sensitivity of c" to scale
lengths, is reported here. The second approach is the subject
of a report by Donaldson [12].
-------
29
Beyond this immediate difficulty, the trial calculations
2
of the second moment of the concentration distribution a
z
exhibit all of the characteristics of classical diffusion
2
theory. As can be seen from the plot of a vs. x in
£i
Figure 4, after an initial period during which the fluxes of
2
material C'w' are brought into equilibrium, o increases
p 2 z 2
as the square of the distance and, as o approaches A ,
this dependence goes over to the classical Fickian diffusion
2 2
rate, with a ~- x . The asymptotic behavior of a is,
Z Z
therefore, analogous to the behavior predicted by the statis-
tical theory of turbulence and we may associate A with
max
the scale length of the turbulence.
If we wish to do so, every term explicit in the model may
be printed out and detailed checks made of their relative
magnitudes and their spatial and time derivatives. To illus-
trate this capacity for detailed analyses of complex situations,
we have plotted the vertical profiles of the turbulent flux of
material at x = 5000, 7000, and 10,000 feet in Figure 5. The
manner and rate at which the initially downward flux below the
plume goes over to an upward flux and establishes a completely
mixed layer is clearly portrayed by these calculations.
With these general features and behavior of the diffusion
model in hand, attention was turned to experimentation with
the only free parameter available, A , and to a first
assessment of the effects of stability and wind shear on the
vertical diffusion predicted by the model.
5.2 An Initial Test of the Effect of Stability
As a first attempt to determine the model's prediction of
the effect of stability on vertical diffusion, the model was
rerun using the same assumptions as before, but inputting the
stable turbulence field as derived from the Hanford data.
Again, the scale length A was set equal to 60 ft, the
iricix
value chosen for the turbulence generation run. The profiles
-------
30
of concentration at x = 500, 5000, and 10,000 feet are shown
in Figure 6 and may be compared directly with Figure 3.
2
The values of a are shown in Figure 7 and, for conven-
2
ience, comparable values of a for the neutral case (Fig.4)
2
are replotted here also. Finally, experimental values of cr
Z
for neutral and stable atmospheres are plotted in Figure 7, as
an initial verification of the model predictions. As can be
seen from these comparisons, surprisingly realistic values of
2
a are obtained for the neutral case. In the stable case, the
model prediction is clearly overpredicting a when compared
2 z
with the Hanford observations, but o vs. x has the right
Z
shape.
It is.interesting to note that the model predicts an
equivalent growth rate for the neutral and stable diffusion
2 2
cases until o reaches approximately 1000 ft . Then the
£1
growth rates diverge rapidly, with the further vertical diffu-
sion of the cloud much slower in the stably stratified
atmosphere. This result suggests that vertical diffusion is
relatively independent of stability and wind shear while the
cloud is small, a result which is observed in the Hanford
experiments and which conforms to the expectation that the
large-scale turbulent eddies are selectively suppressed by
stable temperature stratification (cf., Lumley and Panofsky
[6]).
Considering the sweeping nature of the assumption that a
single scale length and'the mean profiles of wind and temper-
ature correctly model all facets of the structure of boundary
layer turbulence and the diffusion of material embedded in
that turbulence, this initial test can only be considered
remarkable and promising. We can, therefore, undertake further
tests and development of the model with confidence.
-------
31
5.3 The Sensitivity of Diffusion to the Scale Length
Having demonstrated a generally credible simulation
of turbulent diffusion, attention was turned to the anoma-
lous behavior of the diffusion prediction near the core of
the diffusing cloud. For the simplified, two-dimensional
model which we are testing, the equation for the prediction
of the rate of change of local turbulent fluxes of material
(Eq. (14)) may be written
0-C.lw' ~
P^u
p
/^ v"
s-^-/V
P
-(
v/K
0 p
A L
^o +
i
&
'w1
p
^1
p
2
+
s
1
PO
+
^0
a
3
2
•dz
^o V
'^c-
>***
PO
2
p
w1
2
2
Pn
Ai/F P . I
dz y
2C^' ^
xa j
Yap°^sW
(* ) jv
(14)
At the elevation of the core of the cloud (200 feet above
ground), all of the molecular terms are negligible and, for
the cases tested, C'T' is one to two orders of magnitude
less than • C'w' . In noting the behavior of the model, we
may, therefore, focus attention on the first three terms of
the above equation for p u^C'w'/Sx .
M K o p
The first term, -p w'w'(dc /Sz), is a production term
and represents the creation of flux, or.correlation between
C' and w1, out of the mean gradient of C" and the
vertical component of turbulence. We note particularly that
this production'rate depends upon the scale length A only
through the production of vertical velocity fluctuations
w'w1 , a turbulence field which should be independent (under
our present assumptions of a chemically and radiatively inert
diffusing material) of the distribution of C . The local
production of turbulent flux depends, therefore, only on the
-------
32
scale length associated with the production of turbulent
motions and the local mean gradient of the diffusing material.
The second term, 3(
-------
33
turbulent motions are independent of the material distribu-
tions, this approach poses only complications of interpre-
tation, not inconsistencies.
Without prejudging the exact hybrid nature of A,.ff ,
we may systematically solve the diffusion equations for
various choices of this parameter. Intuitively, we presume
the scale length of the turbulence will represent an upper
limit on the. hybrid -^'ff- > but no rigorous arguments have
been developed. Since the primary anomaly which initiated
this line of inquiry was the slow diffusion of the core of
the cloud, we may examine this feature of the model's
dependence upon A,.,,- first.
Figures 8 and 9 show the predicted values of Cn
pmax
at various distances and for various choices of A^'ff ,
using the neutral and stable cases, respectively. These
results show a strong dependence of Cn on A „ and,
^max an i _
in particular, they reveal an axis of minimum values of C
which depends only on ^ f and is independent of stability.
Mathematically, this result simply states that for each
distance of travel, there is a choice of A „ which maxi-
mizes the model's estimate of the turbulent flux of material
out of the core of the plume, and that this choice of A,.,,
increases with increasing,distance of travel.
This result immediately suggests that the hybrid diffusive
scale length A^'fr should be a function of the scale of the
cloud, the one characteristic of the diffusion process which
also increases with distance. Among the candidate scale
lengths of the cloud which might be considered, only one is
independent of the functional form of the distribution of the
cloud concentration, and this is the standard deviation of
that distribution. As another test of the model, we, therefore,
chose to key A,.,,,, on a , and the model was rerun for both
CLJL1 X Z
the neutral and stable cases under the condition A =
= a < A = 60 ft. The resulting predictions for
Z ""~~ IftclX
-------
2 —
a and C are shown in Figures 10 and 11. As can be
seen immediately, this choice for A had only a minor
effect in the neutral case. The decrease of axial concentra-
tion is slightly faster than when A,.,,,, = A^,,^ = 60 feet in
cm i max
the early stages of cloud travel, but under this more energetic
turbulence, diffusion is. still rapid and A,.ff goes to its
limiting value of 60 feet in the first 1000 feet of travel.
A markedly different result is obtained for the stable case,
as is shown in Figure 11. Under this less energetic condition,
cloud growth is markedly suppressed and the axial concentration
is reduced much more quickly during the first 10,000 ft of
2
travel. However, the prediction of a is still too large
£l
beyond 2300 feet, suggesting that A = 60 ft is an over-
max
estimate of this scale length for stably stratified atmospheres.
These results must be considered largely empirical at the
present time and subject to amendment as the results of
further parameter searches become available (cf., Ref. 12).
We may note at this stage, however, that beyond distances of
the order of 10,000 feet, the predictions of the invariant
model, of the classical eddy diffusivity models, and of the
statistical models are indistinguishable for the relatively
simple flow and stability patterns adopted for these tests.
All predict the classical Fickian diffusion of the cloud. At
the shorter distances, however, the invariant model tends to
the statistical theory with only a minimum assumption with
regard to the dependence of the vertical diffusion coefficient
on scale length. No major . assumptions or appeal to empiricism
are necessary in specifying the effective eddy diffusivity as
a function of height and fluid flow conditions. This distinct
advantage is, of course, directly attributable to the retention
of a large measure of the physics and dynamics of turbulent flow
processes in the invariant model. The value of this more
rigorous treatment when applied to arbitrary flow and stability
conditions - conditions for which empirical relationships
become questionable or simply nonexistent - is obvious .
-------
35
6. A CALCULATION OP THE TIME REQUIRED FOR TURBULENCE TO
BECOME FULLY DEVELOPED
A problem of considerable interest in practical air
pollution problems is the rate at which the turbulence
structure, and, therefore, the diffusion rates, adjust to
a major change in the atmospheric flow field. Such changes
might be induced by an abrupt change in surface conditions,
such as flow from land to water, or rapidly changing pressure
gradients with accelerating wind speeds. That portion of the
invariant model used to calculate the turbulence field
provides the necessary information for estimating the rate at
which the atmosphere would respond to such large perturbations.
As an initial experiment to this end, we may postulate the
extreme example of an atmosphere which is initially at rest
and has any prescribed lapse rate of temperature over the
height of interest. At time t = 0, this atmosphere is set in
laminar motion, and a perturbing turbulent component of motion
is superimposed. The model is then run in real time steps
until the components of turbulence and turbulent fluxes come
into equilibrium.
Two examples of the time history of the development of
2 2
turbulence, as measured by the maximum values of uf , w1 ,
p p p
u' + v1 + w1 = K, and u'w1, are shown in Figure 12. The
lower part of Figure 12 is for the case of the unstable to
neutral atmosphere used in the diffusion calculations; the
upper plot is for the case of the Hanford stable atmosphere
used earlier.
As can be seen from these plots, the turbulence structure
in the neutral atmosphere achieved equilibrium in approximately
50-100 seconds. The stable atmosphere required more than ten
times that length of time. This result reflects the strong
damping effect of stability on vertical exchange processes and
the order-of-magnitude larger turbulent energy generated by the
2
neutral atmosphere. (Note the difference in scales of u1 , etc.)
-------
36
Although the conditions of this experiment hardly have a
counterpart in nature, this result does suggest a very real
difference in the adjustment of the atmosphere as it moves
from hot to cold surfaces as compared with movement from
cold to hot. The stabilizing effect of motion of warm air
over a cold surface can be expected to greatly retard the
rate of adjustment of the turbulence.
A more realistic experiment (and one for which necessary
programming was started at the end of this phase of work) is
to simulate the effects of surface heat flux and calculate
the rate of adjustment of turbulence structure which is
initially fully developed for the upwind surface temperature
conditions. It is planned to complete this study along these
more realistic lines in the second phase of work.
-------
37
7. CONCLUSIONS
As a first assembly of the principles and methods of
invariant modeling for the simulation of turbulent diffusion
in the atmospheric boundary layer, the results reported here
must be considered successful beyond reasonable expectation.
What has been demonstrated is that, with only prior know-
ledge of the mean wind and temperature profiles and using a
relationship between the macroscale and the dissipation scale
derived for laboratory scale flow, the structure of turbulence
and the diffusion of matter are simulated well within an order
of magnitude of observed values. In most cases, the verifi-
cation is within a factor of two, a result which has been
achieved for other modeling techniques only after much more
effort and resort to empiricism. Perhaps most importantly at
this stage of development of the invariant model, all of the
classical requirements for asymptotic behavior of the diffu-
sion are met.
In the process of developing these initial tests of the
applicability of invariant modeling to atmospheric diffusion,
we have begun to identify the avenues of improvement of these
simulations. By far the most important of these is the
development of a more rigorous basis for estimating the atmos-
pheric scale lengths and the understanding of these scale
lengths as they affect the balance of production, diffusion,
and dissipation of turbulent fluxes. An extensive reexamina-
tion of the scale lengths associated with laboratory jets has
already been initiated (under Air Force and NASA sponsorship).
The results of this work will provide a logical starting place
for further refinement of scale length choices for atmospheric
boundary layer simulation. The objective of this work will be
to develop equations for the scale length which depend only on
the flow and stability conditions. Since good results have
been obtained under relatively crude assumptions regarding
-------
38
these scale lengths (a result which indicates the primary
physics of turbulence and diffusion is accounted for by the
local gradients of atmospheric motions, temperature, and
pollutant materials), we do not expect a large degree of
sensitivity of the scale lengths to these properties.
Their further consideration does offer the possibility of
improved simulation and, most importantly, confident adapta-
tion of the models to arbitrary flow, stability, and pollutant
source configurations.
In the second phase of this development, we propose the
following activities:
1. Completion of the simulation model programming for
three-dimensional diffusion (appropriate to a continuous point
source);
2. The extension of the turbulence generation model to
include the depth of the planetary boundary layer and the
effects of coriolis forces;
3. The exercising of the model, with improved scale length
estimation techniques for arbitrarily stratified planetary
boundary layer conditions;
4. The description of the invariant models specifically
adapted to atmospheric conditions.
-------
39
REFERENCES
1. Hinze, J.O., 1959: Turbulence, McGraw-Hill, New York,
586 pp.
2. Donaldson, Coleman duP., 1971: "Calculation of Turbulent
Shear Plows for Atmospheric and Vortex Motions," AIAA
Journal 10, 1, pp. 4-12.
3. Donaldson, Coleman duP., 1969: "A Computer Study of an
Analytical Model of Boundary Layer Transition," AIAA
Journal ?, 2, pp. 271-278.
4. Donaldson, Coleman duP . and Sullivan, Roger D., 1970:
"Decay of an Isolated Vortex," A.R.A.P. Tech. Memo. 70-6.
5. Donaldson, Coleman duP . , Sullivan, Roger D., and Rosenbaum,
Harold, 1970: 'Theoretical Study of the Generation of
Atmospheric Clear Air Turbulence," AIAA Journal 10, 2,
pp. 162-170.
6. Lumley, John L. and Panofsky, Hans A., 1964: The Structure
of Atmospheric Turbulence, John Wiley & Sons, New York,
239 PP.
7. Rotta, J., 1951: "Statische Theorie Nichthomogener
Turbulenz," Zeitschrif t ' fur Physik, Vol. 129, p.
8.' Priestley, C.H.B., 1959: Turbulent Transfer in the Lower
Atmosphere, University of Chicago Press, Chicago, 130 pp.
9. Pasquill, P., 1962: Atmospheric Diffusion, D. Van Nostrand
Co., Ltd., London, 297 pp.
10. Pandolfo, Joseph P., Atwater,. Marshall A., and Anderson,
Gerald E., 1971: "Prediction by Numerical Models of
Transport and Diffusion in an Urban Boundary Layer," Center
for the Environment and Man, Hartford.
11. Hilst, Glenn R. and Simpson, C.L., 1958: "Observations of
Vertical Diffusion Rates in Stable Atmospheres," J. Meteor
15, 1, pp. 125-126.
12. Donaldson, Coleman duP . , 1971: "A Progress Report on an
Attempt to Construct an- Invariant Model of Turbulent Shear
Flows," AGARD Conference Proceedings No. 95 on Turbulent
Shear Flows, London, September 1971.
-------
0
0
0
400r
300
Z (ft)
200
20 40 60 80 100 A
.2 .4 .6 .8 1.0 c
..02 .04 .06 .08 O.I -
1 I
' 1 W'W'
1
', I
» I
VTM
\ I
v I
\ \
TURBULENCE RUN
COTE 40
Afurb
Cp(x=0) ^j
A,urb.W>
(ft°F/sec)
too
T (Temperature departure from adiabatic lapse rate, °F)
u (Velocity, ft/sec)
w7^ (ft/sec)2
Figure 1. Input values of mean wind speed u , mean temperature
T , and A used to simulate turbulence structure_ of the lower
atmosphere in an unstable-neutral case (u and T taken from
AFCRL data). Predicted profiles of w'_ and w'T'
and the initial concentration profile C (z) used in later
diffusion calculations are also shown. ^
-------
1.0
T
2.0
—T~
3D
-1—
4.0
—i—
.o w'z,K(ft/sec)2
STABLE ATMOSPHERE
400 -
Figure 2. Same as Figure 1 except for a stably stratified
atmosphere as measured at Hanford. K is the total turbulent
222
kinetic energy, u1 + v1 + w1 __j_
the relatively large value of u'
is small.
and is shown to illustrate
even though wf
-------
350r
300 -
250
200
10,000'
z(ft)
LPD COTE 40
LMAX
o
o
o
o"
D
IO
=60'
o
o
o
m
a
10
150-
100-
o
o
m
a
a.
ID
x(ft)
0
500
5,000
10,000
CPmax
1.000
0.844
0.188
0.041
z(ft)
200
198
211
251
a^(ft2)
1
1,629
20,400
36,440
,3/c3/2
0
-0.15
1.3
1.2
vJ
3.0
19.1
5.8
4.5
-2
10'
jj
10°
Figure 3. Predicted vertical profiles of concentration Cp at 500, 5000, and 10,000 feet
downwind from a crosswind line source located at a height of 200 feet in a neutral
atmosphere. Inset table shows central value and first four moments of C (z) .
ru
-------
CM
*
«*.
^.
CM
I05
10"
10'
10'
EQUILIBRATION
— PERIOD FOR
DIFFUSION MODEL
LPD NEUTRAL-
j i
J L
IOV
10'
10'
10*
10'
DISTANCE FROM SOURCE, x (ft)
Figure 4. Variance of the vertical distribution of concentration
a for the neutral case.
-------
bUU
500
400
300
z(ft)
*y f\n
duu
100
"
-
-
-
r j*^
V ^
- ^"^ \
X^-o. \*
x^SOOO'^ "-- A,
i i i i i i ^
-12 -8 -4 C
» »\ ^
\v\ LPD COTE 40
\\
\ \ SOURCE AT z=200'
\ V
\\x
\ XV
H X °^
\ v Xv
\ ^^
\ V x^«v
* X\ V>k^o
1 \ *»
J \~ ~:>
^flt- ^J3_-O— —ip**f>
*
1 ^x = 10,000'
/
? , , ,iii,.
) 4 8 12 16 20 24
w' (xio4 ft/sec)
Figure 5. Profiles of predicted vertical turbulent flux of matter associated
with concentration profiles shown in Figure 3-
-------
350 r
LPD COTE 40
Figure 6. Vertical profiles of concentration at 500, 5000, and 10,000 feet downwind
from a continuous crosswind source located at 200 feet above the surface in the stably
stratified atmosphere.
-------
icr
10'
icr
CM
*•
««•
*•
CM
10
icr
LPD NEUTRAL
TVA EXPERIMENT
(NEUTRAL)
EQUILIBRATION
— PERIOD FOR
DIFFUSION MODEL
STABLE
HANFORD
EXPERIMENT
(STABLE)
10'
10"
10'
DISTANCE FROM SOURCE, x (ft)
Figure 7- Variance of the vertical concentration distribution
a 2 for both the stable and the neutral cases, and comparison
with experimental measurements.
-------
I02
Amox(ft)
10'
10'
AXIS OF M.N.
'
10
-3
/ /?!!
yio.croo' ,' / ' \\\
! 111
'
U*
/ 5000'
X /' 2500'
i\. «' rS
' \\
a / 1
\ \ N
\
\ \ \ » sv/ / IOQ/
\\N\V /
10
-2
' ' i iN.1 i i 1
10
-I
10'
Cp
max
Figure 8. Summary of LPD predicted values of Cp , the axial
concentration, as a function of constant
(neutral atmosphere).
max
and distance x
-------
Amax(ft)
IU
10'
•"^
- AXIS OF MIN.
&
7 Cpmax. x = IO,000'
\s' ,,500 (y
/ S
"J / 2500'"
.N. / /
r 1 \/ /
: i r\ /
I b A ^V0
x v X /
\ \ l \ /
\ \ v X/
\ \ \ /N
<> A 6 it
\ X X I
\ . \ \ \
i r>° 1 1 1 1 1 1 1 1 1 1 1 V 1 Vi 1 VJ V»
' / 1 *
s / 1 /'
A P VD
xX /"
/ / l\
f .6 J7p
,s / /
' ' /
/ /
1000' /
/ 500'
/ /
/W 4""
/
/
/
\J*
x
i \
,4, . XI i i , , 1
10
-3
10'
10"
Cp
max
Figure 9. Summary of LPD predicted values of Cp , the axial
concentration, as a function of constant A and distance x
(stable atmosphere). max
-------
'19
-------
50
erf (ft2)
Pmax
DISTANCE FROM SOURCE, x (ft)
Figure 11. Comparison of predicted o^ and C^ when
z Pmax
AQV = a < 60 feet for the stable case.
-------
51
100 200 300 400 500 600 700
TIME FROM ONSET OF TURBULENCE GENERATION (SEC)
Figure 12. Time history of invariant model calculations of
22 222
turbulence parameters (maxima of u1 , w1 , K = u' +v' +w' ,
and u'w') showing time required for these parameters to come
to equilibrium in neutral and stable atmospheres.
-------
APPENDIX A
Theoretical Study of the Generation of
Atmospheric Clear Air Turbulence
-------
Reprinted from AIAA JOURNAL, Vol. 10, No. 2, February 1972, pp. 162-170
Copyright, 1972, by the American Institute of Aeronautics and Astronautics, and reprinted by permission of the copyright owner
A Theoretical Study of the Generation of Atmospheric-Clear
Air Turbulence
COLEMAN DUP. DONALDSON,* ROGER D. SULLIVAN,! AND HAROLD ROSENBAUMJ
Aeronautical Research Associates of Princeton Inc., Princeton, N.J.
This paper concerns a numerical study of the time history of the growth and decay of turbulence in the
atmosphere. The study was motivated by a desire to better understand the behaviour of shear-generated turbulence
in the Earth's atmosphere—the problem of clear air turbulence (CAT). The variant modeling technique is used
to derive a closed set of equations that govern this phenomenon. The resulting set of equations allows, for the first
time, a step-by-step numerical calculation of the growth of turbulent disturbances in atmospheric shear flows and
the dependence of this growth on atmospheric instabilities.
Nomenclature §
= heat-flow vector
A.B.C = const
rt2, a3 = const
cc^i'j.Cj.Q = const
cf = specific heat at constant pressure
f = unknown function
a
9iJ
K
k
kr
I
P
= acceleration due to gravity
= the metric tensor
=
= conductivity
= eddy conductivity
= Prandtl's mixing length
= pressure
= K"2 = "2
Presented as Paper 70-55 at the AIAA 8th Aerospace Sciences Meet-
ing, New York, January 19-21, 1970; submitted May 4, 1971; revision
received August 23, 1971. This research was supported by NASA under
Contract NASW-1868.
Index categories; Atmospheric, Space, and Oceanographic Science;
Jets, Wakes, and Viscid-Inviscid Flow Interactions.
* President and Senior Consultant. Associate Fellow AIAA.
t Vice President and Consultant. Member AIAA.
I Associate Consultant; presently Senior Consulting Scientist, Avco
Systems Division, Avco Corporation, Wilmington, Mass. Member
AIAA.
tj The notation of general tensor analysis is used. In particular, the
summation convention is implied, and a comma preceding a subscript
indicates cov.iri.iiu differentiation.
T
t
u, a, w —
Uj =
X =
\,y,: =
AT, =
(5 =
A A =
temperature
time
x, y, and z components of velocity
velocity vector
body force
Cartesian coordinates
general coordinate system
boundary-layer thickness or scale of mean motion
Kronecker delta
",y +»;,,•
function defined in Eq. (66)
scalar measures of length
first and second coefficients of viscosity
eddy viscosity
function defined in Eq. (64)
density
stress tensor
gravitational potential
Subscript
0 = equilibrium conditions of the atmosphere
Superscripts
= departure from equilibrium
= mean
= dep.irtnre Ironi ilic mean
-------
FEBRUARY 1972
GENERATION OF ATMOSPHERIC-CLEAR AIR TURBULENCE
163
Introduction
FOR many years, those engaged in the study of turbulent
motions have made use of the concepts of eddy viscosity and
eddy conductivity to enable them to calculate the conduction of
momentum and heat in these motions. In an incompressible
medium to which these notions have been extensively applied,
the eddy-conductivity concept relates the turbulent flux of
momentum and heat to the gradient of the mean velocity and
temperature fields according to
rij = M"(j + «;,i)-P<"i'u/> = 0* + /'rX«ij + «j.i) (U
and
If we write the velocity vector ut and the pressure as the sum
of a mean part and a fluctuating part whose time average is
zero, i.e., ui = u{+u! and p — p + p' with u.' = p' = 0, Eq. (3) can
be written, after time-averaging, as
S_
Hi
(5)
-q. = /cT,-pcp = (k + kT)Tt (2)
The eddy viscosity HT and the eddy diffusivity kT in these
expressions have, in the past, been determined either by resort
to experiment or, more indirectly, from physical arguments which
determined some basic and perhaps simpler parameters which
could in turn be evaluated experimentally. The latter method
has, in general, been preferred because of the hope that such
parameters might have some simple general relationships by
which they could be determined for a large class of flows. Of
course, both these methods fail if iy or qi are not zero when
HJJ+ Uj, and T( vanish.
In view of the fact that the tensor and the vector
both have dynamics of their own quite apart from
(although related to) the dynamics of the mean fields u; and T,
the concept of eddy transport has been remarkably successful
over the years in dealing with so many engineering problems.
The explanation for this success lies in the fact that the dynamics
of the turbulence is fast enough relative to the dynamics of
the mean flow in many cases so that the turbulent correlations
and <«/7"> are able to "track" the development of the
mean motion. This is certainly true of developed pipe flows
where the derivatives of all mean quantities, except the static
pressure, with respect to streamwise positon are zero.
On the other hand, a number of turbulent flows of interest
exist in nature for which the turbulent correlations are not able
to track the changes in the mean field uf and T. In such flows,
the dynamics of the turbulence is important, and the concept of
eddy transport in any simple form is not a sufficiently powerful
tool with which to study the motion.
As will be seen in what follows, the motion of the atmosphere
is an example of such a flow. For this reason, a more powerful
approach must be used—one which actually keeps track in some
way of the dynamics of the turbulent correlations such as
and .
Within the past few years, a method has been developed by the
senior author of this paper and his collaborators1'2 which will,
when it has been fully developed, permit one to calculate, with
some degree of precision, the motion of strongly sheared
turbulent layers. The method, which is called the method of
invariant modeling, has been applied with some success to the
calculation of incompressible turbulent boundary layers and to
the decay of an isolated vortex. In this paper, we present the
results of an initial application of the method to the calculation
of atmospheric motions.
Brief Description of the Method of Invariant Modeling
It is beyond the scope of this paper to present a complete
description of the method of invariant modeling. However, we
give, in this section, a brief outline of the method, as well as a
short discussion of its relationship to the mixing length model
of turbulent transport proposed by Prandtl.3
Consider the motion of an incompressible medium which is
governed by the Navier-Stokes equations so that the momentum
equation may be written
d(pu^ldt + (puiu,)J = - P,j + i\j (3)
where
t;: = n(u, j+u,() (4)
This is the turbulent momentum equation derived by Reynolds.
To solve this equation we need to know the Reynolds stress
tensor — p. If the turbulent correlations can track the
mean motion in the sense discussed in the previous section,
we can assume that this tensor can be expressed in terms of
the mean velocity field u,. •
Prandtl accomplished this for the case of a simple boundary-
layer flow by means of a physical argument3 and obtained the
well-known result
- p = p/2 |rti/<>| nl /'<> (6)
where / (the "mixing length") is a scale parameter identified in
the argument with a characteristic size of the turbulent motion.
We can obtain this result in another way. We ask whether one
can write the tensor in terms of the mean motion u,.. The
expression chosen must exhibit the symmetric-tensor character
of , must be dimensionally correct, must be invariant
under a general transformation of coordinates, as well as under
a Galilean transformation, and must not be such as to upset any
of the general conservation laws. Since
eii=«i.i+"i,i <7)
has the required symmetry and satisfies the invariance conditions,
one might write
-<«/«/> =/(«*K"ij+»j.i) <8*
where /(ut) must be a scalar and must be such that the
dimensions of the resulting expression are correct. Two simple
choices for f(uk) have the required Galilean invariance. They
are
./>*) = A2", mm (9)
and
In these expressions, A is a scalar measure of length. The first
expression above is zero by virtue of continuity, so we choose
the second expression [Eq. (10)] so that
- = A2(elnBfi'"")l/2(iJu+ «,.,.) (II)
If the formula given above is reduced to the situation of a simple
shear layer u = u(y) with v = w = 0, one obtains
- = (2)1/2A2 \du/i)y\ du/Py (12)
which is essentially Prandtl's result. We conclude, therefore, that
Eq. (11) represents a general form of Prandtl's mixing length
expression for the Reynolds stress. We might expect such a result
to hold provided the mean flow was such that the dynamics of
the turbulence could track the mean motion. We also expect the
scale parameter A to be related to the scale of the mean motion.
Suppose now we wish to generalize this idea to the case where
the dynamics of the turbulence must be considered. We can
write the equation for the Reynolds stress, namely,
P<""""/"/>.m- <";Y>,j- <«/P'>,i+ +
^mn<<";>,™-2/' (13)
and we note that correlations other than the second-order cor-
relation itself appear in the equation.
If we wish to close the set of equations to be solved with no
more equations than those describing u. and , we must
express these extra terms through the use of quantities for which
we already have equations. In the method of invariant modeling,
we do this by choosing expressions for these terms based on the
second-order correlations themselves which satisfy the general
rules for modeling set forth above in the discussion of Prandtl's
-------
164
DONALDSON, SULLIVAN, AND ROSKNBAUM
A1AA JOURNAL
mixing length. When this procedure is followed, one can actually
obtain only a few models which satisfy the rules set down and
which are reasonably simple expressions representing straight-
forward physical principles. The following model has been
investigated in some detail:
/> .] (14)
(15)
X'>.m (17)
In this model there are two scale parameters A and A. They are
related in concept to the integral and microscales of the turbu-
lence. In our work, we have taken them to be related to each
other by the expression
A2 = A2/(Cl + c2«c'A) (18)
Equation (18) has been used by Glushko4 and by Beckwith and
Bushnell5 when modeling the turbulent kinetic energy equation
obtained from Eq. (13) by contraction.
In applying the model to the incompressible boundary layer
on a flat plate, we assume that in the outer regions of the
boundary layer the scale A is a constant fraction of the
boundary-layer thickness 6
A = c3«5 (19)
Near the surface in the boundary layer we assume that A is pro-
portional to the distance normal to the surface, i.e.,
A = c4y (20)
for v -» 0. An examination of the boundary conditions at y = 0
shows, however, that the constant c4 must be related to f, so
that only three parameters, f ,, c2, and c3, must be determined
from experimental data.
To see the relationship between c, and c4, we note that at the
surface we must have
meters that are presently being generated are inserted into our
atmospheric programs.
Equations of Atmospheric Motion
Since the variation of density with altitude plays an important
part in atmospheric motion, the equations above do not apply.
The appropriate equations are generally written7 in terms of the
departure of the various physical quantities, 7* p. p, and u., from
the values they would have in an atmosphere at rest, namely,
T0, p0, p0, (which are functions of altitude only) and u, = 0.
Thus, we define
f = T- T0; p = p- p0; p = p-p0; and ui = t<( (29)
These quantities, in turn, are split into their mean and fluctuat-
ing parts, so that, for the rest of this paper, T stands for the mean
of f and T = t— t. The other variables are treated the same
way.
With the assumption of small departures from equilibrium, the
equations for the turbulent atmosphere are
T)/rf + (f<0ujf)j =
(30)
(31>
(32)
(H -f /(*) [(Po .m/Po)j + <""""/> (Po .m/Po) .,- +
(Po,>o)« ",'"j"">
(34)
Now must be zero at y = 0 and, since /«? through the surface y = 0 by viscous
diffusion which is not physically possible, we must have
r<«('M/>//> = «2>'2 + a3y3-K.. (22)
Substituting Eq. (22) into Eq. (21) results in
A2 = 2(a2 y2 + a3 y3 + . . .) /(2a2 + 6a3 y) (23)
or (unless u2 — 0) for y -> 0
A = y (24)
Since Re -> 0. Eq. (18) can be written
A2 = A2/f ,
for small y. Thus, from Eqs. (20) and (24) we get
(25)
y2 = c4V/c,
or
(26)
(27)
We have applied this model to the turbulent boundary layer
on a flat plate and have compared the results with the standard
data on such flows6 in an attempt to determine the best values
off,, c2, and c'3. When this was originally done, the following
values for the parameters were found
= 0.125; c3 = 0.064
(28)
Subsequently, an error was found in the computer program so
that the values of the parameters given previously are not now
considered to be those that will give an optimum fit to existing
experimental data. They are, however, not far different from the
optimum values. In the work that will be presented in the
remainder of this paper, the values given in Eq. (28) were used.
We believe, therefore, that the results we shall present are correct
as to the general behavior of the atmosphere, although the exact
numerical values of the various quantities computed will be
altered somewhat when the new "optimum" values of the para-
(35)
(36)
cjk.
that
In these equations, we have taken the Prandtl number, n
equal to one, but we have not made the usual assumption
u J' = 0 as is done in Ref. 7 and which, in general, restricts the
motions which may be studied to those which are small in scale
compared to the scale of the atmosphere.
If Eqs. (34, 35, and 36) are studied, one sees that, if the
terms containing p0~ ' p0 f are dropped, the invariant modeling we
have already carried out for boundary-layer flows may be taken
over without the introduction of new parameters to the calcula-
tion of atmospheric motions. If we keep the terms containing
Po~ 'PO ,•• we nnd tnat tne modeling is changed slightly so as to
satisfy the conservation law given in Eq. (33), The important
point, however, is that no new parameters are introduced and the
three parameters f ,, c2, and c3 discussed in the previous section
are sufficient to determine the motion. We use the following
modeling
<"/w;>.n + <«""i',-'>J+.i] (37)
Ml'r>ill+j (38)
= -\qcr= -
-------
FEBRUARY 1972
GENERATION OF ATMOSPHERIC-CLEAR AIR TURBULENCE
165
Wj +u'jj}> = (p0]-
«P'> Po,/Po - <"/P'> Po,i/Po <45>
= -(PO«/A)<«,T'> (46»
«« /"> = - <«""»/> p0>o (47)
=-<«"'r >/;0>0 (48)
Equations (30-36), together with the modelings given above, now
constitute the equations with which we propose to study the
nature of turbulence generation in the atmosphere as a result
of winds and thermal instabilities. Of course, this can only be
done for very simple motions for which we will be able to
identify the scale of the mean fields of u( and T In what follows,
we will apply the equations given in this section to a particularly
simple example of atmospheric motion with the idea of demon-
strating some of the behavior of atmospheric shear layers as
indicated by our equations.
A Simple Case of Atmospheric Shear
In order to study the nature of atmospheric motion as
described by the equations given above, we have selected what
we consider to be the simplest motion that can give meaningful
results. We consider an atmosphere in which the motion (u, v, w)
in Cartesian coordinates (\. y, :) is such that u = u (:, t) and
f = u- = 0. Thus, we assume that the mean motion consists of
only one component of velocity. That component of velocity is
parallel to the ground and is only a function of time and
altitude. In such a motion, symmetry requires that the correla-
tions , . and be zero. To gain this amount of
simplicity in the motion, it is necessary to invent a body force
X(:.r) which starts the motion. Actually, there is no such real
force acting on the atmosphere but rather the atmosphere is
driven by pressure and Coriolis forces. These forces, however,
produce a far more complicated motion in the large than the one
we have chosen. At the very least, one must, in such cases, con-
sider a motion of the form u(z, t) and w(z.t). Nevertheless, these
real motions can resemble locally the simple motions we shall
study here. We believe that the essential features of atmospheric-
turbulence production will be demonstrated by our somewhat
fictitious atmospheric model while retaining a relatively simple
set of governing equations. For the motion we have just
described, the equations given in the previous section become
P0?u/dt = (/ifVf'z2)-['W"'vO)/'--] + A'(r,t) (49)
Po df/ct = (n P2f/dz2)- [rlp0O'r»A\-] (50)
rp/rr= -(ryfzXp0) + (p00/70)f (51)
; , fiu 1 f
_.= _2 — + — T
f
ft
/
A(
i fV
, 1 p
Po<:\
<«v> *)
3/
riA
0 -
3
- ; - = — — IP
ft Po<
Po
. ,
— - 1 V>_-
\
(52)
,53)
-.2
Po '"-*
Po '• Po
LU ff/J+T0
Po '• Po
(55)
= _---+- -^U ,,A)-
x f'z fz pn fz V (": J
<\w'T">
- __
ft
Po'-
2/i
Po *
\
vw> — + — —
(z p0ez
(56)
Po
(57)
(58)
Po''" \Po '-
From our previous work on invariant modeling, we will keep the
following results
/.2 = A2/(2.5 + 0.125 ReA) (59)
A = 0.064^ (60)
In this last equation, we will assume that f> is some appro-
priately chosen scale of the mean motion which will be akin to
the boundary-layer thickness of our previous studies.
In the present study, we do not consider any interaction of the
turbulence that is produced with the ground since we will, in
general, be studying the generation of shear layers at altitude.
However, if we were to do so, we would choose
A=(2.5)"2.- (61)
in the region from the ground to the point where
_- = [0.064/<2.5)"2]<} (62)
Initial Computations
To date we have made only a few computations bused on the
scheme laid out in the preceding section. These initial computa-
tions represent the simplified cases which were used to check the
program and are not meant to represent any particular physical
situation. As may be seen by an examination of Eqs. (49 57),
if we neglect the terms containing = \r'2> =
= 1 (fps)2. The band is centered at an altitude of 20,000 ft
although the actual value of the altitude has no bearing on the
problem due to the assumption that p0 is constant. (The initial
sharp edges of the turbulent band are soon smoothed off by the
action of the diffusion and dissipation terms of the equations.|
The atmosphere is initially at rest, i.e., at t = 0, u = 0. For a time.
a body force acts on the atmosphere, creating a mean motion.
The body force in the calculations we present here was taken
to be of the form (Fig. 1)
X(:.() = <•(I -c2)2 (63)
-------
166
DONALDSON, SULLIVAN, AND ROSENBAUM
AIAA JOURNAL
Fig. 1 Shape of the mean-
velocity forcing function.
where
for
= (r-20,000)/1000
(64)
<
(65)
In Eq. (64), z is an actual altitude in feet. Outside the region
|£| < 1, the driving force is taken to be zero. Thus, the forcing
function is such that it produces a mean shear layer some 2000 ft
thick. In the particular examples which we consider here, this
forcing function is such that in the absence of any turbulence or
viscosity the centerline velocity would increase to 30 fps in
3000 sec. After 3000 sec, it was assumed to vanish and the
atmosphere left to dissipate the motion by turbulence and viscous
diffusion.
Four cases of initial mean temperature distribution were con-
sidered (Fig. 2).
Case I: neutrally stable atmosphere T(z,0) = 0.
Case II: unstable-stable-unstable atmosphere given by
f(z,0) = (1 35/32)^(2 -|£|)2 for |c| < 2 and zero elsewhere. Notice
that this initial temperature is spread over 4000 ft rather than
2000 ft as was the case for the forcing function. The maximum
temperature difference is 10°R.
Case III: stable-unstable-stable atmosphere given by
T(z,0)= -(135/32)£(2-|£|)2 in the region]^ < 2 and zero else-
where.
Case IV: stable-unstable atmosphere given by T(z,0) =
\0(\-ri2)2 where
r, = if = (z - 20,000)/2000 (66)
for \t]\ < 1 and zero elsewhere. This formula again gives a
maximum temperature difference of 10°R.
The other initial conditions that have been used for these
computations are
<«V>(z,0) = 0; (2,0) = 0; and (z,0) = 0
The scale A for all calculations was taken to be the length
defined by 0.064 times the breadth of the mean shear layer as
based on the distance between the points where the mean
velocity falls to half its maximum value.
Altitude, kilofeet
22
CASE HI
Fig. 2 Initial temperature
profiles.
Fig. 3 u. , and K = + + <>'2> for a neutrally stable
atmosphere; Case I at 300 sec.
Results
Case I : Neutrally stable atmosphere
In Figs. 3, 4, and 5, we show the calculated values of the
mean velocity u, K = <«'2> + + , and for three
times, 300, 3000, and 6000 sec, respectively. The first figure shows
the very start of the development of the mean motion. At this
time, very little turbulence has been created by the mean motion.
Most of the turbulent energy that is seen is left over from the
decaying initial turbulent layer. In Fig. 4, which shows the motion
at the end of the application of the forcing function, the mean
velocity is a maximum. The production of turbulent energy as a
result of the deformation cu/?z is evident. The magnitudes of
K and that are calculated are consistent with the results of
tests on freejets in the laboratory. Figure 5 shows the motion
after decay has set in. The mean velocity has spread, and the
centerline value of K has started to fill in by diffusion.
Another picture of the development of this motion can be had
by plotting the maxima of the various correlations determined
at any time against time. This is done in Fig. 6. Here we see the
initial decay of the turbulent layer with which we started, fol-
lowed by the production of the various correlations as a result of
the interaction between this turbulence and the mean flow. Once
the forcing function ceases to increase the mean velocity, the jet
that has been formed decays in the typical self-similar manner
of laboratory jets.
Case II: Unstable-stable-unstable atmosphere
In Fig. 7, we show the mean-velocity profile at 3000 sec for
the production of a shear layer in an unstable-stable-unstable
,f!2/sec2
Fig. 4 a, , and K for Case I at 3000 sec.
-------
FEBRUARY 1972
GENERATION OF ATMOSPHF:RIC-CLEAR AIR TURBULENCE
167
,f12/sec2
-2 -I
K,ft2/sec2
•(u'w'^.fl /sec
2, 2
Fig. 5 u, . and K for Case I at 6000 sec.
Fig. 7 it and for Case II at 3000 sec.
atmosphere. We note from a comparison of this figure with Fig. 4
(and the result is generally true) that at the maximum of the
mean velocity, the mean velocity and shear — are
changed somewhat but not drastically by the effect of the tem-
perature profile. As may be seen in Fig. 8, where we have
plotted K and for both 300 and 3000 sec, the develop-
ment of the turbulence is significantly altered. At small times, the
contribution of the vertical component of velocity to the
turbulent energy is large and is due to the basic thermal
instability. The maximum value of K at this time is well outside
the region of normal production by tu/Pz. As time progresses,
however, the production of turbulent energy by f u/r.: becomes
large, and the maximum of K occurs near the maximum of
ru/rr as may be seen from the K curve for 3000 sec.
In Fig. 9, we show the behavior of the heat-flux correlation
and the mean temperature for times of 300 and 3000 sec.
Here we note the very large values of in regions of
cf/d: <0 and the very small values of in regions of
fT/f: > 0. The values of |rT/rr are not very different, and the
level of turbulence K at times greater than 3000 sec cannot
account for this difference. That is to say, none of the usual
eddy diffusivity models could account for this behavior and the
very persistent nature of the inversion that has been calculated
here. A close inspection of the solution in the neighborhood of
z = 20.000 (the center of the inversion) indicates that at small
times (t less than approximately 300 sec), the solution for
10°
t, seconds
In Fig. 10, we have again plotted the behavior of the maxima
of the velocity correlations as a function of time. Note first the
general production of turbulence by thermal instability at small
times. Near the maximum of the mean velocity, the production
of turbulence resembles that of the neutral atmosphere. Once the
forcing function disappears, the mean motion dissipates and the
thermal instabilities in the outer regions of the motion take over.
In this case they just about balance dissipation as time proceeds.
This is due to the increase in the dissipation scale /. in our
calculations as time proceeds.
Case III: Stable-unstable-stahle atmosphere
In Fig. 11, we show the mean velocity at 3000 sec for the
production of a shear layer in a stable -unstable -stable atmo-
sphere. Note again that this is not much changed from the t\\o
previous cases although there has been slightly more diffusion of
the mean velocity profile.
In Fig. 12, we again see a marked effect of the mean-
temperature profile on the turbulence production. In this case, the
thermal instability at small times produces a great deal of turbu-
lence in the unstable layer trapped between two stable layers.
Only towards the maximum of the mean motion does the shear-
produced turbulence make its presence fell.
Figure 13 is a plot of the behavior of and 7" Here we
see again that very large turbulent heat fluxes are developed in
the unstabe region cT/f: < 0, while it has been virtually impos-
sible to produce any heat flux in the region where ff/i: > 0.
This is another example of the difficulty of getting rid of inver-
300 seconds
3000 seconds
Fig. 8 <7"2> and K
for Case II at 3000
sec.
Fig. 6 Maxima of velocity correlations; Case I.
-------
168
DONALDSON, SULLIVAN, AND ROSENBAUM
AIAA JOURNAL
Alti
kilo
24
23
22
21
10 -5
i i _J
1 Jt
i.o ° Vy
\\ 19
ude,
feet
3000 seconds
I
^X 5 10 15
&^**^ j I 1
0.5 1.0 1.5
\ ,ft°R/sec
"V
300 seconds
— 3000 seconds
\
1 '
17'/
i
16 -
Fig. 9 f and <»v'T'> for Case II at 300 and 3000 sec.
10
Velocity
correlations, 1-0
f(2/sec2
O.I
t, seconds
Fig. 10 Maxima of velocity correlations; Case II.
18
I7L
Fig. 12 < r 2> and K for Case III at 300 and 3000 sec.
sions. Again this reluctance of nature to establish a heat flux
when f-T/c: > 0 is a basic characteristic of the atmosphere and
cannot be described by any simple eddy-difTusivity model.
In Fig. 14 we show the behavior of the maxima of t lie-
velocity correlations as a function of time for this case of shear-
layer formation. Note that the effect of shear-generated turbu-
lence makes itself evident only when the mean flow is near its
maximum velocity.
Case IV: Stable-unstable atmosphere
In Fig. 15, we have plotted the mean-velocity profile and the
profile of for a time of 3000 sec. Here we again note
the similarity of this profile to all the previous cases.
In Fig. 16, we show the behavior of < T2> and K for
r = 3000 sec. In this case we note that there are two extrema
of K. Both are close to the point where ru/fr is a maximum, as
might be expected. The production of turbulent energy on the
unstable side exceeds that on the stable side but not by as
much as we had expected prior to making these calculations.
In contrast to the behavior of the turbulent energy, we see
from Fig. 17 that the production of in a stable region, so
that the temperature inversion (rT/rr > 0) is extremely persis-
tent. Again, no presently used eddy-diffusivity model would
predict the result we have shown.
In Fig. 18, we show the behavior of the maxima of the
velocity correlations as a functon of time. We note again the
initial production of turbulence by means of thermal instability,
the modification of this by the action of shear-generated turbu-
lence near the time of maximum mean velocity, and the final
phase of the motion when the production of turbulence by
300 seconds
w'T'>, ft«R/sec
Fig. II u and for Case III at 3000 sec.
Fig. 13 f and for Case III at 300 and 3000 sec.
-------
FEBRUARY 1972
GENERATION OF ATMOSPHERIC-CLEAR AIR TURBULENCE
169
Velocity
correlations, i.o
Fig, 14 Maxima of velocity correlations; Case III.
Fig. 16 <7'2> and K for Case IV at 3000 sec.
thermal instability in the decaying thermal layer is just about
balanced by dissipation so that a rather constant level of
turbulent energy results.
Discussion of Results
The most interesting aspect of the results we have just presented
is the radically different behavior of the heat flux correlation
depending on whether dT/dz is greater or less than zero.
The explanation for this behavior lies in the equations and may
be exhibited as follows.
Consider a flow in which there is no mean motion and in
which all quantities except Tare independent of z. In this case,
the pertinent equations for the production of correlations
obtained from Eqs. (54, 57, and 58) are
Pt
Tn
Po
dt
d (68)
(69)
If we consider- only the production terms and neglect the ten-
dency-towards-isotropy terms and dissipation terms in these
equations, we have
a/3t = (2g/T0) (70)
d\)
(72)
,ftz/sec2
-a -i
300 seconds
3000 seconds
1.5 2.0 2.5
,ft°R/sec
Fig. 17 f and for Case IV at 300 and 3000 sec.
10
Velocity
correlations, i.o
ft2/sec*
O.I
Xv'V>
10 3
t, seconds
Fig. 15 u and for Case IV at 3000 sec.
Fig. 18 Maxima of velocity correlations; Case IV.
-------
170
DONALDSON, SULLIVAN, AND ROSENBAUM
AIAA JOURNAL
We now assume (on the basis of experience with existing
calculations) that the rate of change of dt/P: is very slow com-
pared to the rate of change of the second-order correlation. In
this case we can assume, to first order, that df/ and from Eq. (71) by differentiating
that equation with respect to time and using Eqs. (70) and (72).
The result is
?2/(]r2 = -(4&/T0)(^fA\-) (73)
Examination of this equation shows that when ft/tlz < 0, there
is an exponential development of . However, when the
atmosphere is stable, i.e., fTff; > 0, the heat-flux correlation
is oscillatory with frequency
« = 2[(0/70)a77rz]1'2 (74)
In fact the general solution for Eqs. (70-72) may be written
= Asinwt + Bcoswt (75)
= [cy/(2^f/r?z)][-/lcoscot + Bsinwr] + C (76)
= (2/to)(cf/?z)[Acos'T>,
, and . In fact, if A is sufficiently large compared
to /. so that only dissipation need be considered, the solution
becomes exp[ — 2/ir/(p0 A2)] times each of the expressions given
above.
A close examination of all our numerical results shows this
oscillatory behavior at the frequency given above during the
early phases of each run in the region where cf/<~>z > 0.
The upshot of all this discussion is that the basic equations
of atmospheric motion indicate a vast difference between the
behavior of when cf/dz < 0 and when ff/f; > 0. In the
unstable case, a large heat flux is established quite rapidly. In the
stable case, the heat flux tends to oscillate about a value that is
very small. This results in an extraordinary persistence for tem-
perature inversions in the atmosphere. It also means that the
idea of a simple eddy diffusivity in the atmosphere is completely
wrong and will, when used in calculations, result in far too rapid
a dissipation of any inversions.
Conclusions and Recommendations
We have presented the first results using the method of
invariant modeling to calculate the behavior of atmospheric
shear layers. We realize we have only scratched the surface; much
remains to be done. In particular, we are working towards a
better method for putting the scale length A into our equations.
We would eventually like to have an explicit equation for A which
could be carried along as an extra simultaneous equation to be
solved. For certain cases, namely, the production of thermal
turbulence, this is probably more important than in the case of
shear-driven turbulence. In this connection, mention should be
made of the recent work of Harlow et al.8'9 which is closely
related to our own studies.
Although the authors recognize that much still needs to be
done to finalize a general method for predicting wind-driven
atmospheric motion, we feel that the method of invariant model-
ing represents a powerful basic tool with which to approach this
problem. It seems clear from our studies to date that the eddy-
diffusivity concept should be discarded as a means by which the
gas dynamicist attacks the problem of turbulent motions in the
atmosphere.
In regard to future applications of the method we have just
presented, it should be mentioned that equations similar to those
considered here can be written for correlations of fluctuations in
species concentration. Thus, application of this method to the
convection and diffusion of gaseous pollutants in the atmosphere
is straightforward. In particular, the effect of a stable lapse rate on
the development of turbulent mass transport is similar to that
found for the turbulent transport of heat in this study.
References
' Donaldson, C. duP., "A Computer Study of an Analytical Model
of Boundary Layer Transition," AIAA Journal, Vol. 7, No. 2, Feb. 1969,
pp. 271-278.
2 Donaldson, C. du P. and Rosenbaum, H., "Calculation of Turbulent
Shear Flows Through Closure of the Reynolds Equations by Invariant
Modeling," Rept. 127. Dec. 1968, Aeronautical Research Associates of
Princeton Inc., Princeton, N.J.
3 Prandtl, L., "The Mechanics of Viscous Fluids," Aerodynamic
Theory, Vol. Ill, edited by W. F. Durand, Dover Publications. New
York, 1963.
4 Glushko, G. S., "Turbulent Boundary Layer on a Flat Plate in an
Incompressible Fluid," Bulk-tin of the Academy of Sciences of the
USSR. Mechanical Series, No. 4, 1965.
5 Beckwith, I. E. and Bushnell, D. M., "Detailed Description and
Results of a Method for Computing Mean and Fluctuating Quantities
in Turbulent Boundary Layers," TN D-4815, 1968, NASA.
6 Coles, D. E. and Hirst, E. A., editors, Proceedings Computation
of Turbulent Boundary Layers, AFOSR-lFP-Stanford Conference,
Vol. II, Aug. 18-25, 1968.
7 Lumley, J. L. and Panofsky, H. A., The Structure of Atmospheric
Turbulence, Interscience Publishers, New York, 1964.
8 Harlow, F. H. and Nakayama, P. I., "Turbulence Transport
Equations," The Physics of Fluids, Vol. 10, No. 11, 1967.
9 Harlow, F. H. and Romero, N. C., "Turbulence Distortion in a
Nonuniform Tunnel," Rept. LA-4247, 1969, Los Alamos Scientific Lab,
Los Alamos, N.Mex.
------- |