Statistical Modeling and Simulation cf
Atmospheric Turbulence
By
SURESH B. PAHWA
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA IN PARTIAL
FULFULLMENT OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1975
ACKNOWLEDGEMENTS
The author wishes to express his sincere appreciation to
Prof. R.W. Fahien, chairman of the supervisory committee, for his
enthusiastic interest and valuable guidance throughout this research.
He is grateful to Dr. D.W. Kirmse, cochairman of the supervisory
committee, for his interest and many valuable suggestions. He wishes
to thank the other members of his supervisory committee, Prof. R.D.
Walker, Dr. K.E. Gubbins, Dr. A.E.S. Green and Dr. W.C. Huber, for
their cooperation in serving on the committee.
The author wishes to thank Dr. A.W. Westerburg for the use
of the least squares subroutine RMINSQ.
The author gratefully acknowledges the financial assistance
of the College of Engineering, University of Florida. He wishes to
express his gratitude and appreciation to the Department of Chemical
Engineering, University of Florida for both financial support and
personal assistance from its faculty and staff.
Finally, the author extends his thanks to his wife, Brigitte,
whose active support made this work worthwhile.
TABLE OF CONTENTS
Page
ACKNOWLEDGEMENTS................... ..... .............. ii
NOMENCLATURE.......................................... vi
ABSTRACT............... ............. ............... xi
CHAPTER:
1. INTRODUCTION................................... 1
2. PREVIOUS STUDIES ............................... 4
2.1 Structure of Atmospheric Turbulence....... 4
2.1.1 Definition........................... 4
2.1.2 Theories of turbulence................ 5
2.1.3 Turbulence measurements.............. 9
2.1.4 The MoninObukhov similarity
theory............................... 12
212 Atmospheric Diffusion...................... 14
2.2.1 Early Developments and the
Ktheory............................. 14
2.2.2 Gaussian plume models................. 18
2.2.3 Some recent studies.................. 19
2.2.4 Urban air pollution models........... 21
2.3 Stochastic Models of Turbulence............ 23
3. THEORY............... ...... ........... ....... 27
3.1 Theory of Stochastic Processes............. 27
3.1.1 The FokkerPlanck equation........... 27
3.1.2 Stochastic differential equation..... 33
3.1.3 Comparison of the FokkerPlanck
equation with the diffusion equation. 37
TABLE OF CONTENTS (continued)
3.2 The Growth Law ...........................
3,2.1 The Lagrangian similarity theory
for large times......................
32.2 The growth law for a Markov process...
THE STOCHASTIC MODEL............................
4,1 Description of the Model.....................
4.
4.2.2.5
4.2.2.6
The mean time steps...........
The mean Lagrangian
accelerations.................
. RESULTS..........................................
5.1 Computational Parameters....................
Estimation of Simulation Parameters.........
4,2.1 Neutral conditions....................
4.2.1.1 The mean velocity profile.....
4.2.1.2 The variance of the velocities
4.2.1.3 The variance of the
accelerations.................
4.2.1.4 The mean time steps...........
4.2.1.5 The mean Lagrangian
accelerations.................
4,2,2 Nonneutral conditions................
4.2.2.1 The mean velocity profile.....
4.2.2.2 The temperature profile.......
4.2.2.3 The variance of the velocities
4.2.2.4 The variance of the
accelerations................
4,2
Page
40
40
44
52
52
55
59
59
61
TABLE OF CONTENTS (continued)
Page
5.2 Validation of the Model..................... 85
5.2.1 Neutral conditions.................... 86
5.2.2 Nonneutral conditions ................ 96
5.3 Correlated Results.......................... 107
5.4 Concentration Distribution in the
Vertical Direction.......................... 144
5.4.1 Additive function for deviation from
Gaussian behavior...................... 149
5.4.2 NonGaussian exponential function..... 153
6. CONCLUSIONS AND RECOMMENDATIONS.................. 159
6.1 Conclusions... ............................. 159
6.2 Recommendations... ......................... 161
APPENDICES:
A. DISCRETIZATION OF DISTRIBUTION CURVES AND
THE MONTE CARLO SIMULATION METHOD............... 163
B. COMPUTER PROGRAMS.................................... 170
LITERATURE CITED....................................... 206
BIOGRAPHICAL SKETCH...................................... 222
NOMENCLATURE
2
a acceleration, Lt2
a constant in Eq. (4.52)
b constant in Lagrangian similarity theory (Sec. 3.2)
B ratio of the Eulerian and Lagrangian integral length
scales
c concentration, ML3
c concentration distribution coefficient, ML3
3
c ground level concentration at y=0, ML
go
2 2 1
C specific heat at constant pressure, L t T
2 1
D molecular mass diffusivity, L t
E expected value (Sec. 3.1)
E energy spectrum function,L t2
E Gaussian concentration distribution function in the
y ydirection
f deterministic function (Sec. 3.1)
f space correlation coefficient (Sec. 4.2.2.4)
f probability density function (App. A)
F characteristic function of turbulence (Sec. 2.1.4)
2 1
F energy spectrum function, L t
F concentration in the zdirection
z
g acceleration due to gravity, Lt2
g random function (Sec. 3.1)
h deviation function from Gaussian distribution in the
zdirection without reflection
h' reflection term corresponding to h
H stack height, L
H
z
I
k
K
L
L
L
c
LG
ml, mi2
M
n
n
n
N1
N2
Ncell
N
P
p
p
Pl' P2
P
q
q1' q2
Q
 deviation function from Gaussian distribution in the
zdirection
 concentration integral, ML2
 wave number, L1
2 1
eddy diffusivity, L T
 length scale of turbulence, L
 MoninObukhov length scale, L
 dimension of length
 length scale for plume spread, L
gradient similarity length scale, L
 exponents in Equations (3.45) and (3.46)
 dimension of mass
1
 frequency, t
 partition number for cumulative distribution function
(App. A)
 exponent in Sutton's Lagrangian correlation coeffici
ent expression (Sec. 2.2.1)
 stability parameter E TLG/6R
 dimensionless pressure gradient
 number of particle samplings in a cell
 total number of particles released
 transition probability
 exponents in Equation (3.40) and (3.41)
1 2
 pressure, ML t2
heat flux, Mt3
 constants in Equations (2.14) and (2.15)
source strength, Mt1.
 source strength, Mt
r generation of the species due to chemical reaction,
ML3t1
r exponentially distributed random variable with mean
of one
R gas constant (Sec. 4.2.2.6), MI2t2T1
R correlation coefficient
s normally distributed random variable with mean=0, and
variance=l
t time
tE Eulerian microscale of turbulence, t
t sampling time, t
T temperature
TG gradient similarity temperature scale, T
T MoninObukhov temperature scale, T
T integral time scale of turbulence, t
l
u velocity in the xdirection, Lt1
S friction velocity, Lt1
Up convection velocity, Lt
v velocity in the ydirection, Lt
V velocity vector, Lt1
V c volume of a cell L
w velocity in the zdirection, Lt
x coordinate axis, downwind distance, L
X displacement of a particle in the xdirection, L
X displacement vector, L
y coordinate axis, crosswind distance, L
Y displacement of a particle in the ydirection, L
viii
z coordinate axis, vertical distance, L
z roughness length, L
Z displacement of a particle in the zdirection, L
Greek Letters:
a velocity profile exponent
a nt moment of particle displacement, Lntl
S ratio of the Lagrangian and Eulerian integral time
scales
S constant in the loglinear velocity profile (Sec. 4.2.2)
8 exponent in Eq. (5.60)
Y exponent in Eq. (5.56)
r adiabatic lapse rate, TL1
6 exponent in F equation (Sec. 5.4)
z
2 3
e rate of turbulent energy dissipation, L t
n independent variable in the Lagrangian similarity
theory = dt/T
6 potential temperature, T
+ dimensionless potential temperature at the ground
eR
K von Karman constant
A mean value of the time steps
t
A integral length scale of turbulence, L
S viscosity, ML tl
2 1
v kinematic viscosity, L t
dimensionless vertical distance (Sec. 5.4.2)
p density, ML3
a standard deviation
T lag time, t
T time scale of turbulence (Sec. 3.2), t
1 2
S shear stress (Sec. 4.2.2.1), ML t
1 2
T shear stress at the ground, ML t
dimensionless velocity scale
w velocity scale of particle motion, Lt
Subscripts:
a acceleration
E Eulerian
h heat transfer
i i direction
L Lagrangian
m momentum transfer
o initial values
o steady state neutral conditions
R reference
t time
v velocity
S' wind direction
Superscripts:
fluctuating quantity
mean quantity
* dimensionless variable
deviation from steady state neutral conditions
t turbulent
Abstract of Dissertation Presented to the Graduate Council
of the University of Florida in Partial Fulfillment of the Requirements
for the Degree of Doctor of Philosophy
STATISTICAL MODELING AND SIMULATION OF
ATMOSPHERIC TURBULENCE
By
Suresh B. Pahwa
March, 1975
Chairman: R.W. Fahien
Cochairman: D.W. Kirmse
Major Department: Chemical Engineering
A statistical model of turbulent diffusion as applied to flow
in the atmosphere is presented in this work. The model was originally
developed and applied to flow in a pipe by D.W. Kirmse and R.W. Fahien
in 1964.
It is shown that the displacement of a particle as modeled in
this work can be represented as a stochastic process with Markovian
properties. The differential equation for the velocity of the particle
is compared with a generalized stochastic differential equation for a
Markov process. Since the displacement of a particle is assumed to be
Markovian, it obeys a FokkerPlanck equation. Comparison of the Fokker
Planck equation with the diffusion equation showed that the deterministic
velocity used in the FokkerPlanck equation in general is not equal to
the mean wind velocity, and that the Lagrangian diffusivity is equal to
the Eulerian diffusivity only under certain circumstances.
An approach similar to G.K. Batchelor's similarity theory
was used with the assumption of a Markov process to arrive at a
growth law that expresses the spread of a plume in the vertical
and crosswind directions, as it moves downwind.
In the simulation, a large number of fluid particles were
released from the source and followed downwind as they move according
to the convection and diffusion in the atmosphere. The particle
displacements are sampled at a fixed interval of time in the downwind
cells to obtain the pollutant concentrations. The particle displace
ments were simulated by integrating the turbulent Lagrangian accelera
tions. The total acceleration is the sum of a mean and a fluctuating
acceleration. Estimates of the mean Lagrangian accelerations were
calculated from the equations of motion in the Lagrangian form with
the molecular effects neglected. The initial fluctuating velocities,
the fluctuating accelerations and the time steps were obtained from
the desired probability distribution curves. The statistical parameters
used in the distribution curves were estimated from the turbulent in
tensities, the mean velocity, the ratio of the Lagrangian to Eulerian
integral time scales, and either the turbulent energy dissipation or
the auto correlation curves.
A sensitivity analysis was carried out to test the effect of
the internal parameters of the model, the total number of particles
released, the sampling time, and the cell size. Optimal values of
these parameters were used in the runs.
The mean Lagrangian accelerations for neutral conditions were
taken as zero. The results obtained from the present work are com
pared with the experimental atmospheric data of Project Prairie Grass
and with predictions made by other methods.
xii
For nonneutral conditions, the concentration distributions
were found to change with variations in the value of the mean Lagran
gian acceleration in the vertical direction. However, reasonable
agreement with the experimental data is observed when the mean Lagran
gian accelerations are taken as zero.
The correlated turbulence data of Cramer and Record for neutral
conditions were used to study the effects of the stack height, the
roughness length and the friction velocity. The results were then cor
related to obtain generalized charts in dimensionless form for the eddy
diffusivities and the dispersion coefficients as functions of the down
wind distance with the dimensionless roughness length as the parameter.
xiii
CHAPTER 1
INTRODUCTION
The achievements of the industrial world during the present
century have been accompanied by great ecological disturbances which
have affected the quality of life on the earth. The air we breathe
today is more contaminated than what it was even a decade ago. Eco
logically balanced urban designs require that the harmful effects of
air pollution be minimal. In the event of a nuclear fallout, one
should be able to predict the spreading of the radioactive material
with sufficient accuracy and speed. Since the Clean Air Act of 1967,
some serious efforts have been made in the abatement of industrial
pollution.
The goals of reaching certain pollution levels are related
to the design of emission control equipment through the use of a dif
fusion model. A good diffusion model should be able to correlate
pollutant concentration as a function of the rate of emission, source
height, the ground roughness, and meteorological conditions. A sta
tistical model for turbulent diffusion in the atmosphere from a
steady state, continuous, point source is presented in this work.
Turbulence is stochastic in its behavior. It is not possible
to predict the detailed structure of turbulence at any given time.
Statistical tools enable one to describe turbulence in a comprehen
sive manner. The stochastic model of turbulence as employed in this
work is capable of using the measured or correlated statistical
properties of turbulence and predicting the diffusion in the atmos
phere with reasonable accuracy.
The model was originally developed in a Ph.D. thesis by
Kirmse (1964) with Fahien, and applied to turbulent flow in a pipe.
Later Gielow (1972) used the model for flow between parallel plates.
In this work, the model is applied to the atmospheric ground layer.
The model uses a Lagrangian point of view. A Monte Carlo method is
used to generate random values of the Lagrangian accelerations which
are integrated twice to get the fluid particle displacements. A
large number of particles are released from the source and their
paths simulated until they leave the region of interest. The pollu
tant concentration is obtained by sampling the particle positions in
the cells down wind of the source.
As is shown in Section (3.1), the displacement of the particle
is modeled as a stochastic process with Markovian properties, and thus
it obeys the FokkerPlanck equation which can be considered as a form
of the diffusion equation using the Lagrangian point of view. The
correlation coefficient for a Markov process is of the exponential de
cay type, and according to Taylor's (1921) theory, an expression for
the standard deviation of the concentration distribution can be derived.
The model was verified with the field observations for differ
ent stability conditions and reasonably good agreement was observed
with the data. The correlated turbulence data were used to improve
the PasquillGifford charts for the neutral conditions so as to include
3
the effects of the source height and the roughness length. Also,
the eddy diffusivity components were calculated; these can be used
in the solution of the diffusion equation.
CHAPTER 2
PREVIOUS STUDIES
2.1 Structure of Atmospheric Turbulence
2.1.1 Definition
In general, fluids are subjected to two kinds of random
movements: (i) thermal agitation of molecules and (ii) turbulence.
The molecular random movements are discontinuous due to the collisions
of molecules with each other, and with the wall. The molecules have
well defined shapes and sizes which are not dependent upon these
movements. In contrast, turbulence is a continuum motion, and the
eddies change their shapes and sizes continuously. Sutton (1949)
compares the contrast between laminar flow and turbulent flow as
being similar to the contrast between a musical note and a noise in
acoustics.
It is difficult to define turbulence. Perhaps the most
complete definition available in literature is given by Hinze (1959,
page 1): "Turbulent fluid motion is an irregular condition of flow
in which the various quantities show a random variation with time
and space coordinates, so that statistically distinct average values
can be discerned." The researchers in general have found it easier
to describe turbulence than give a formal definition for it.
Bradshaw (1971) describes turbulence as a tangle of vortex lines or
partly rolled vortex sheets in a preferred direction by the mean
flow and in random directions with each other. Corrsin's (1961)
article is a good introduction to turbulence. He describes turbulence
as a random field of vorticity, and states that turbulence can be
expected in a fluid whenever there is a shearing flow and the inertial
effects are much larger than the viscous effects. Lumley and Panofsky
(1964) opted to list the essential features of turbulence. Their
description can be summed up as: Turbulence is a mode of fluid motion
that is a threedimensional continuum phenomenon, rotational, dissi
pative, nonlinear, stochastic and diffusive with the transport
occurring at time and length scales typically the same as those of
the properties being transferred.
Other introductory discussions of turbulence are given by
Sutton (1953), Schlichting (1955), Townsend (1956), Brodkey (1967),
Batchelor (1967), Phillips (1969), Monin and Yaglom (1971), Tennekes
and Lumley (1972), and Vinnichenko et al. (1973).
2.1.2 Theories of turbulence
Historically, Reynolds (1883) was the first person to visualize
turbulence in a laboratory. His classical experiment consisted of
injecting a dye in a glass tube and observing its mixing behavior.
The smooth flow of the dye with the mean flow changed to rapid lateral
mixing when a critical Reynolds number was reached.
Prior to 1921, turbulence was thought of as consisting of
discrete motions and collisions in an analogous manner to the
molecular motion (Batchelor, 1967). G. I. Taylor (1921) broke with
the idea of analogy to the kinetic theory of gases and introduced
the correlation between the velocities at two points as one of the
quantities needed to describe turbulence. His theory also clearly
shows the enormous simplifications obtained by assuming homogeneity
and isotropy. This assumption has been widely used since to achieve
simplifications in formulation of the turbulence theories.
Richardson (1922) suggested that the energy from the mean
flow is drawn by the large eddies and the large eddies in turn break
up into smaller eddies until the smallest size of eddies, i.e., the
dissipation scale is reached.
Keller and Friedmann in 1924 were the first to state that the
correlation functions and the statistical moments are the fundamental
characteristics of turbulence (Monin and Yaglom, 1971). Using the
equation of motion of a real fluid, Friedmann proposed a general
method of obtaining the differential equations for the statistical
moments of arbitrary order. The total system of equations consists
of an infinite number of equations. Any finite subsystem is always
nonclosed. Turbulence modeling techniques associated with the
closure problem are of major interest to investigators in the field
of turbulence.
Based on an analogy to the kinetic theory of gases, Prandtl
(1925) proposed a mixing length of turbulence similar to the mean
free path of the molecules. The theory effectively states that the
turbulent kinematic viscosity is the product of a local characteristic
velocity and the mixing length.
Vt a v'Z
(2.1)
Taylor (1938a) showed the existence of interactions between
components having different length scales. Taylor (1938b) studied the
energy spectrum function of turbulence. The energy spectrum function
describes the distribution of turbulent kinetic energy over various
eddy sizes. He showed that this energy spectrum function is the
Fourier transform of the correlation between the eddy velocities at
two points. Von Karman (1938a and 1938b) showed that the velocity
correlations are tensors. This has enabled more precise and concise
working equations for the nonisotropic situations. Von Karman and
Howarth (1938) reduced these tensors to scalar quantities in the case
of isotropic turbulence. This exemplifies the simplification achieved
through the assumption of isotropy.
Kolmogorov in 1941 put forward the theory of universal
equilibrium. It states that the turbulence energy is generated in
a certain range of length scale and then passed on to smaller and
smaller scales by a cascade process. Thus the sizes are only in
directly coupled. At high Reynolds numbers, the creation and the
dissipation ranges are widely separated (Ellison, 1962). Obukhov
(1941) arrived at the same conclusions independently. Kolmogorov
arrives at the concept of local isotropy. As the sizes of the
eddies get smaller, there is a tendency for the turbulence structure
to be isotropic. Kolmogorov's theory leads to the conclusion
(Tennekes and Lumley, 1972) that the turbulence energy in the overlap
region of the energy creation range and the Kolmogorov range can be
expressed as
E(k) a 82/3 k5/3 (2.2)
where E is the energy spectrum function, E is the rate of turbulence
energy dissipation, and k is the wave number. This range is known as
the inertial subrange.
Kolmogorov (1962) and Obukhov (1962) refined their theory of
the small scale structure depending upon one value of the energy
dissipation E (and v) only. In fact, E is a fluctuating quantity,
and the mean E varies with the Reynolds number. A slight deviation
from Eq. (2.2) is observed by taking into account the variations in
the energy dissipation.
Kraichnan (1959) proposed a solution to the closure problem
by expanding the third and higher order statistical moments of the
Fourier transform of the velocity in terms of the second order moments
and their response functions. The response function of the time
correlation coefficient with the lag time T is exp (vk2 ), where k
is the wave number. He applied the theory to stationary, isotropic
turbulence of very high Reynolds number, where the Reynolds number
was defined in terms of the wave number and the turbulent intensity.
The energy spectrum function in this range was obtained in terms of
a nonlinear integral involving the time correlation function and
response functions.
Various closure theories have been put forward by Heisenberg
(1948), Kovasznay (1948), Proudman and Reid (1954), Chandrasekhar
(1955), Tatsumi (1957), Patankar and Spalding (1967), and others.
A survey of developments in the field up to 1970 is discussed by
Reynolds (1974).
2.1.3 Turbulence measurements
One of the more significant sets of experimental measurements
was made by Laufer (1954) for air in fully developed flow in a pipe
at Reynolds numbers from about 50,000 to 500,000. The mean flow and
the statistical properties measurements include the mean velocity
profile, radial, axial and angular turbulence intensities, cross
correlations between the different components of the fluctuating
velocity, energy spectra and the turbulent kinetic energy dissipa
tion rate.
In the atmosphere, Richardson (1919,1921,1926) did not use
any of the refined instruments available today, but was able to make
some meaningful approximations and arrive at useful results. Taylor
(1935,1936) studied the behavior of isotropic turbulence in a decaying
wind stream.
Taylor (1938b), using a hot wire anemometer, measured the
kinetic energy associated with different wave numbers, i.e., the
energy spectra. Deissler (1950) studied adiabatic turbulent flows
in pipes. Favre et al. (1953,1957,1958) measured space and time
correlation functions in a wind tunnel and a boundary layer.
Lauffer (1950) conducted a detailed study of the mean and fluctuating
quantities in a two dimensional channel. He made the energy balances
and verified Kolmogorov's hypothesis.
The isotropic turbulence was the object of the studies of
Liepmann and Lauffer (1951), Hanratty (1956) and Corrsin (1958).
In a study of turbulence in a boundary layer over a flat plate
Klebanoff and Diehl (1952) and Klebanoff (1955) found that the
probability densities of the fluctuating velocities are very close
to Gaussian. Other important contributions in the measurements of
laboratory scale turbulence include the works of Rotta (1950),
Liepmann and Robinson (1953), Pai (1953), Sandborn (1955), Sandborn
and Braun (1956), Hanratty (1958), Frenkiel (1958), Johnk and Hanratty
(1962a,1962b), Barbin and Jones (1963), Cermak and Baldwin (1963),
Baldwin and Mickelsen (1963), Clark (1968), Grass (1971), and Snyder
and Lumley (1971).
Sheppard (1946) measured the velocity profile in the atmosphere.
Sheppard used a drag plate at the ground to measure the drag force.
The velocity profile was observed to be logarithmic for near neutral
conditions.
1 f rz+zo
= ( log ] zz (2.3)
where K is the von Karman constant, T is the shear stress on the
ground, p is the density of air, and zo is the roughness length.
The value of the drag coefficient as a function of the roughness
length was evaluated and the variation of the mixing length with
height was studied.
For one observation site, Tait (1949) obtained an empirical
relationship for the temperature gradient in terms of the intensity
of radiation. Lettau (1959) examined the effect of surface roughness
on the geostrophic drag coefficient. The geostrophic drag coefficient
is the ratio of wind speed to the geostrophic speedwhere the coriolis
force balances the pressure gradient force. He also discussed the
effect of Richardson number on the shear stress and the velocity
profile. Panofsky and Deland (1959) analyzed the onedimensional
energy spectra in the lowest one hundred meters of the atmosphere.
They found the Taylor hypothesis regarding equivalence of Eulerian
time and space correlations to be correct, and observed that the
Lagrangian and Eulerian correlations have similar shapes. Record
and Cramer (1966) compared the energy dissipation rates obtained
by previous investigators and found that the variations in c and with
stability and height were extremely large. The turbulence intensity
and the energy spectra data measured at Round Hill were correlated
in the later studies (Cramer et al., 1967; Cramer and Record, 1969).
Other measurements of the turbulent intensities and the energy spectra
include the works of Taylor (1952), Gifford (1955,1956), Panofsky et al.
(1958), Thompson (1962) and Angell (1964).
The study of fluxgradient relations and the temperature
fluctuations in the ground layer have been conducted by Lettau (1950),
Corrsin (1951), Bolgiano (1959,1962), Priestley (1960), Lumley (1964),
Webster (1964), Dyer (1965), Gurvich and Zubkovskii (1966), Clarke
(1970), Okamato and Webb (1970), and Haugen et al. (1971). Haugen
et al. (1971) experimentally found that momentum and heat fluxes are
constant to within 20% provided that a sufficiently long averaging
period is chosen. The velocityvelocity and velocity temperature
correlations were measured and it was shown that the vertical velocity
fluctuation w' is better correlated with temperature than with the
horizontal velocity fluctuation u'.
L
The simulation of the atmospheric ground layer in a wind tunnel
has been of some interest to the experimentalists in recent times.
Strom and Halitsky (1955) outline various factors to be considered in
designing a wind tunnel. Lloyd (1967) generated a shearing flow in
a wind tunnel with an artificially roughened floor. Pruppacher and
Neilburger (1968) discuss the design of a wind tunnel at the University
of California. Arya (1968) and Schon et al. (1974) similated thermally
stratified boundary layers in wind tunnels and compared their observa
tions with the atmospheric data. Cermak et al. (1966) performed a
dimensional analysis for a wind tunnel and numerous experiments.
2.1.4 The MoninObukhov similarity theory
The similarity theory of Monin and Obukhov (1954) by dimensional
analysis arrives at characteristic length and temperatures in a ther
mally stratified gound layer.
Vertical temperature gradients in the atmosphere create density
differences large enough for the buoyancy term to be included in the
equation of motion. However, based on empirical observations, the
vertical heat and momentum fluxes can be assumed to be constant in
the ground layer. The other assumptions that must be made to derive
the similarity theory are: (i) complete homogeneity exists in the
horizontal (xy) plane, (ii) variations of density produced by pressure
fluctuations are neglected, (iii) the temperature deviations from the
neutral conditions are small and (iv) the viscous effects are neglected.
Then the three components of the equation of motion, the equation of
continuity, and the thermal energy equation reduce to:
Du 1 8P
Dtu p Px (2.4)
Dt p ax
Dv 1 8P
Dv al (2.5)
Dt p ay
Dw 1 p
Dt p + g T (2.6)
Dt p z TR
Bu av aw
S+ a + 0 (2.7)
ax ay az
where u, v and w are the velocities in the x, y and z directions
respectively, p is the density of air, P and T are the pressure and
temperature deviations respectively from the neutral conditions, g is
the acceleration due to gravity, and TR is the reference temperature.
The basic hypothesis is that the above equations along with
the assumption of the constancy of the heat and momentum fluxes
describe the statistical characteristics of the turbulence in the
ground layer completely. The characteristic parameters for the layer
are then: g/TR, u* and q/C p where u, is the friction velocity. The
only length and temperature scales that can be formed by a combination
of these parameters are:
3
U*
L = (2.8)
T C P
T = (2.9)
Ku, Cp
If the similarity theory holds, the dependence of any mean
turbulence characteristic F on height may be written in the form
F = i (z/L) (2.10)
o
where F is a parameter with the dimensions of F, formed by combination
of g/To, p, u, and q/C p, and t is a universal function of z/L.
Panofsky's (1961) analysis based on the measurements of
Priestley (1960) shows that the similarity theory holds for the
velocity and temperature profiles. The similarity theory was used
by Charnock (1967) to correlate Richardson number, dimensionless
velocity and temperature gradients as functions of the dimensionless
height z/L. The Richardson number can be interpreted as the gravita
tional force over the inertial force. Taylor (1960), Klug (1968) and
Swinbank (1968) compared the observations made in the atmosphere with
the similarity theory.
Calder (1968) is critical of the similarity theory, and says
that the similarity theory is not a strict consequence of dynamical
equations. There is no closed system of equations for any finite
number of statistical characteristics of turbulencenot even for the
velocity and the temperature gradients. The hypothesis of the existence
of only three parameters, g/TR, u*, q/C p, is a specific viewpoint on
the mechanism of turbulence, and is a hypothesis of maximum simplicity.
2.2 Atmospheric Diffusion
2.2.1 Early developments and the Ktheory
The time averaged equation of continuity for a species in a
15
medium can be expressed in general as (Bird et al., 1960):
Dc 2
+ VV'c' = DV c + r (2.11)
Dt
where c is the time averaged concentration of the species, V' is the
fluctuating velocity vector, D is the molecular diffusivity, and r is
the rate of production of the species by chemical reaction. A simple
form of Eq. (2.11) is obtained by making the following assumptions:
(i) the turbulent mass flux can be expressed in a manner
analogous to that used by Boussinesq (1877) to describe the turbulent
momentum flux; i.e., it is taken as being proportional to the potential
gradient or
V'c' = K*Vc (2.12)
(ii) molecular diffusion is neglected,
(iii) chemical reaction is neglected.
Roberts (1923) used Eq. (2.12) assuming the eddy diffusivity
to be a scalar quantity (which implies isotropy) that is independent
of position, for the cases of dispersion of a pollutant from (i) an
instantaneous point source, (ii) a continuous point source in a uniform
wind, and (iii) a continuous infinite line source perpendicular to the
direction of a uniform wind. He also obtained an expression for the
visible outline of a smoke plume, and extended the results obtained
with an isotropic diffusivity to anisotropy.
Richardson (1926) considered the diffusion process as the
scatter of a cluster of pollutant particles. He used a scalar dif
fusivity that was independent of the position in space but was a
function of the effective eddy size by considering the distance
between the diffusing and the neighboring particles. Sutton (1932)
developed the idea of the effective eddy size. He started with
Taylor's (1921) theory
t t'
2 2 f f
2 = 2c f J RL(T) dT dt' (2.13)
0 0
where a2 is the variance of the dispersion, 2 is the variance of the
v
velocity, and RL is the Lagrangian correlation coefficient. He assumed
that
n
RL() = 2 (2.14)
V +a T
n > 0
and derived an expression for a. His result showed that for n < 1, the
size of the effective eddy increases with the distance of the diffusing
cluster from its source.
Sutton (1934) developed an expression for the mass flux due to
evaporation from a surface assuming that the diffusion of mass follows
the same laws as the diffusion of momentum. Sutton (1947) extended his
previous theory of eddy diffusion to include the turbulence intensities
and extended the analysis to heat transfer as well (Sutton, 1948).
Bosanquet and Pearson (1936) proposed a solution of Eq. (2.11)
which is well known as the Ktheory. Velocity is assumed independent
of height and the eddy diffusivities in the y and the zdirections
are assumed as follows:
Ky = q u x (2.14)
K = q2 u z (2.15)
where K and K are the eddy mass diffusivity components in the y
y z
and zdirections respectively, q1 and q2 are constants, u is the
constant mean wind velocity, and x and z are the mean wind and the
vertical distances respectively.
Calder (1949) approximated the logarithmic velocity profile
for neutral conditions by the power law profile, and obtained the solu
tion of the twodimensional diffusion equation. Calder (1952) made a
contribution to the Ktheory by expressing the eddy diffusivity as
K = K u, z (2.16)
z
where u* is the friction velocity, and K is the von Karman constant.
Starting with a power law form of the velocity profile and
assuming the turbulent Schmidt number to be one, Deacon (1949) arrives
at the following form for the vertical turbulent diffusivity:
Kz = K u* z (2.17)
a > 1 unstable conditions
a < 1 stable conditions
Smith (1957) also used a power law form for the vertical eddy
diffusivity and assumed that the eddy diffusivity in the y
direction also varies in an identical manner with height. He
obtained expressions for ground level concentration profiles from
elevated sources. Further developments of the Ktheory include
the works of Davies (1950,1954), Knighting (1952) and Rounds (1955).
2.2.2 Gaussian plume models
Hay and Pasquill (1957) measured dispersion down wind of a
source and showed that the Lagrangian correlations are maintained
high for much longer periods of time than the Eulerian correlations.
This work formed the basis, and later led to the development of the
Gaussian plume models. Hay and Pasquill (1959) explored the idea
further and hypothesized that the Lagrangian and the Eulerian correla
tions are similar in shapes but have different time scales. Implying
a correlation coefficient equal to one, they assumed
cy Og x (2.18)
where a is the standard deviation of the concentration distribution
y
in the ydirection (cross wind direction), Oa is the standard deviation
of the horizontal wind direction and x is the down wind distance.
However, the correlation coefficient in practice is always less than
one, and Eq. (2.18) should be modified accordingly. Based on the
available data (Haugen, 1959; Monin, 1959; Smith and Hay, 1961; and
a few others), Pasquill (1961) plotted ay vs. x for different 0
values. Similar curves for 0z vs. x were also plotted. In the same
z
paper, Pasquill also showed how to estimate Oa from qualitative
meteorological parameters.
Gifford (1957) showed how similarity predictions of relative
diffusion theory (Batchelor, 1949,1952) apply to the spreading of
smoke puffs, and, using the observations of Kellog (1956) and Frenkiel
and Katz (1956), correlated mean square particle dispersion as a func
tion of time. Gifford (1959,1961) estimated the diffusion coefficients
in Sutton's (1932) formulation and outlined the method to estimate the
parameters needed to use Pasquill's (1961) graphs. Pasquill and
Gifford both worked independently to estimate the a's and the charts
originally formulated by Pasquill are better known as PasquillGifford
charts. Bowne (1973) has presumably improved the charts by including
the effect of the down wind terrain in a qualitative manner.
The Gaussian plumes are widely used in industry and by the
pollution control agencies (Frenkiel, 1956a; Davidson, 1967; Croke,
1968; Shieh, 1969).
2.2.3 Some recent studies
The number and types of the variables involved in an atmospheric
dispersion problem make the problem a very complex one in nature.
Panofsky (1969) summarizes the different meteorological and physical
parameters involved and their effects. Seinfeld (1969) includes the
effect of chemical reactions and surveys the different diffusion models.
Introductory books on the subject are those of Sutton (1953),
Priestley (1959), Pasquill (1962), Slade (1968), Csanady (1973), and
Seinfeld (1975). The emphasis in the current research work has
shifted to the prediction of dispersion rather than presentation of
theories. This has led to a certain amount of field measurements
and to the estimation of eddy diffusivities.
Stewart et al. (1958) used radioactive argon as a tracer
material from a 61 m stack (BEPO reactor) to study the diffusion of
the plume in the atmosphere down wind of the stack. The effective
values of the constants in Sutton's (1932) diffusivity expressions
were evaluated. Saffman (1962) studied the effect of wind shear on
horizontal spread from an instantaneous ground source. Islitzer and
Dumbauld (1963) showed that the deposition of particulate pollutants
can be significant, thus reducing concentrations at greater distances
to well below that predicted by diffusion models. Fuquay et al.
(1964) (Green Glow program) carried out 46 ground source experiments
and concluded that the standard deviation of the wind direction aa
is not a parameter to use for cross wind diffusion, but the product
u0 should be used instead.
Blackadar (1962) related the turbulent momentum diffusivity
K to wind shear through the use of the expression
m
K = 1/3 4/3 (2.19)
m
where k is a length scale of turbulence. Starting from the turbulent
energy equation, the energy dissipation E was expressed in terms of the
wind shear for neutral conditions. Wu (1965) extended the work for
nonneutral conditions, and, assuming the thermal diffusivity Kh to
be proportional to K derived an expression for Kh.
Hanna (1968) hypothesized that the vertical eddy diffusivity
K depends upon the parameters that determine the energy spectrum
of the vertical velocity fluctuations. The observations of Record
and Cramer (1966), Kaimal (1966), Smith (1967) and others were used
to calculate the momentum, mass and thermal diffusivities. Hosler
(1969) estimated the vertical diffusivities from measurements of
atmospheric radon concentrations.
The modeling of chemical reactions and'incorporating them
into the diffusion equation is discussed by Lamb (1973), Heines and
Peters (1973a,1973b) and Fabrick (1974).
The source height used in the dispersion models is the
effective stack height which is the sum of the physical stack
height and the plume rise. Numerous empirical plume rise models
are available in literature, among which Holland's equation (Thomas
et al., 1963) is the most widely used. Moses and Storm (1961)
compared observed plume rises with some wellknown formulas.
Briggs(1965) from dimensional analysis arrived at different equa
tions for transitional, stable and neutral conditions.
2.2.3 Urban air pollution models
The pioneering work in the field was done by Frenkiel (1956b)
who developed a model for the city of Los Angeles. Lucas (1958) in
troduced the use of the Gaussian plumes. This approach was adopted
in several modeling studies: Pooler (1961) studied the nature of
atmospheric diffusion over St. Louis, Clark (1964) modeled the city
of Cincinnati, Turner (1964) developed a model for Nashville, and
Miller and Holzworth (1967) for New York City. Koogler et al. (1967)
included a first order reaction rate in Miller and Holzworth's model,
and applied the model to the city of Jacksonville.
Since then various approaches and numerical techniques have
been employed: Panofsky and Prasad (1967), Davidson (1967), Slade
(1967), Sandberg et al. (1970), Eschenroeder and Martinez (1970),
Takamatsu et al. (1970), Johnson et al. (1971), Wayne et al. (1971),
Roth et al. (1971), Lantz et al. (1971), Bowne et al. (1971) and
Aziz et al. (1973). Comprehensive summaries and discussions can be
found in Stern (1970) and Fan and Horie (1971).
The classifications of Moses (1969) and Roth et al. (1971)
are used here to categorize the urban air pollution models.
(i) The source oriented models: The models of this kind use
the integrated form of the diffusion equation and require the
dispersion coefficients as functions of distances from the source,
a complete source inventory and information about other physical
and meteorological parameters.
(ii) The receptor oriented models: The contributions to pollutant
concentrations at a monitoring station are modeled. The sources may
be autos, industry, space heating or power plant stacks. These models
are useful when a large number of sources of different nature are
involved. Several types of receptor oriented models have been developed.
Among these are the Clark model (Clark, 1964), the regression model
(Roberts et al., 1969) and the Argonne tabulation prediction scheme
(Moses, 1969).
(iii) The diffusion equation models: These models are based on
the numerical solution of the diffusion equation including the
chemical reaction term. These models are further divided into:
a) moving cell approach
b) fixed cell approach
In the moving cell approach, concentration changes in a
hypothetical parcel of air are calculated as it is convected with
the wind. The wind velocity is assumed constant with height, and
the diffusion in the horizontal plane is neglected.
In the fixed cell approach, the airshed is divided into a
threedimensional grid, and the numerical solution may use a finite
difference scheme, a variational method or a particleincell (PIC)
method. The PIC method involves moving the mass particles according
to convection and diffusion.
2.3 Stochastic Models of Turbulence
One of the first stochastic models was developed in a Ph.D.
thesis by Kirmse (1964) and presented by Kirmse and Fahien (1967).
The probability distributions of the random forces acting on the
fluid particles in a random flow were calculated. If the mass of a
fluid particle remains constant, then according to Newton's second
law, the force is proportional to the acceleration. Therefore, it
is sufficient to calculate the probability distributions of the
particle accelerations. The particle's movements are described by
the Lagrangian accelerations as it moves with the flow. The experi
mental data of Liepmann and Robinson (1953) indicates that the
turbulent velocity and its time derivatives have approximately normal
distributions. This information is used to calculate the statistical
parameters and the paths of the fluid particles are simulated on a
digital computer. This model is used in the present work and is
described in detail in Section 3.1.
Kirmse (1964) simulated the turbulent field in a pipe using
the model, and calculated concentration profiles and eddy diffusivity
components downstream of a source. Gielow (1972) used the same model
for flow between two parallel plates and showed that the model worked
for heat transfer as well. Gielow (1972) also simulated the velocity
record and correlation curves and thus showed that the model was
capable of predicting the detailed structure.
Kirmse's (1964) work must be recognized as a pioneering work
in the numerical simulation of an extremely complicated flow structure
of turbulence. It should also be noted that the KirmseFahien model
is not an exact representation of a turbulent flow field. The
Lagrangian velocity and acceleration distributions are assumed to be
of the same nature. A fluid particle is assumed to maintain its
identity throughout, and the entrapment of mass by the eddies is
neglected. It is also assumed that an equal amount of mass of the
diffusing species is carried by all the eddies. It is however quite
obvious from the results that these effects are small enough to be
neglected in studying the phenomenon of turbulent diffusion.
Davies et al. (1954) and Gupta (1959) noted that Taylor's
(1921) representation of particle diffusion by discontinuous movements
implies zero correlation between the velocities during two different
time steps, and is therefore an inadequate representation of turbulence
as was also indicated by Taylor (1921).
Kraichnan (1970) simulated an isotropic velocity field in a
turbulent flow. The velocity field was stored in a computer as a set
of Fourier components, which in turn were picked from Gaussian distribu
tions. Lagrangian velocity correlation, eddy diffusivity and response
functions were calculated and compared with the results obtained by
the direct interaction approximation method (Roberts, 1961).
Deardorff and Peskin's (1970) numerical model involved obtaining the
eddy velocities at each time step from the equation of motion, and
these were integrated to obtain the paths of the particles. The
Lagrangian correlation and mean square dispersion were calculated.
Bullin and Dukler (1974; Bullin, 1972) used a numerical
model similar to the Brownian motion model for the molecular trans
port, and used a white noise signal to represent the molecular motions.
A major difference between the KirmseFahien and the BullinDukler
stochastic models is the representation of the eddy velocity. An
eddy velocity field continuous in shape is obtained from the
KirmseFahien model (Gielow, 1972) and it possesses the desired
properties. While Bullin and Dukler carefully take into account the
molecular action, they represent the eddy velocity field by a rather
crude approximation of a Gaussian white noise. A white noise consists
entirely of uncorrelated impulses adjacent to each other (Jenkins and
Watts, 1969). Its auto correlation function is a dirac delta function
which has a zero auto correlation at all times t / 0 and is infinite
at zero time lapse. The impulses are randomly distributed along the
time axis according to a Poisson distribution (Seinfeld and Lapidus,
1974). The model is similar to the Brownian motion model and may
represent molecular transport adequately, but is an oversimplification
of the complicated nature of turbulence.
Goldstein (1951) presented a different analysis of diffusion
by discontinuous movements, although he considered equal time steps
and a velocity of constant magnitude. His model is similar to
Taylor's (1921) analysis of diffusion by discontinuous movements,
with the difference being that Goldstein introduced a two time step
correlation affecting the sign of the velocity. This assumption
results in a telegraph equation in contrast to a first order
differential equation for the one step time dependent process.
However, it is possible to have a nonzero correlation for the one
step dependent process if the velocity and the time steps are not
of constant order of magnitude. Bourret (1960,1961) uses a similar
method to that of Goldstein (1951) but considers the threedimensional
movement of a particle. Thus he arrives at a sixth order partial
differential equation.
In general, one may expect an equation of infinite order if
the movement of the particle depends on all its previous steps. Such
an analysis is of course inconvenient to deal with. Even the most
complex behaviors in nature, such as the birth and death process and
the genetic effects, can be adequately described by a Markov process
that leads to first order differential equation known as FokkerPlanck
equation.
CHAPTER 3
THEORY
3.1 Theory of Stochastic Processes
A stochastic process is a collection of random variables, or
a chain of events that follow a certain probability law (Tucker,
1967; Kahn, 1956). The random variables can be considered as the
values of a single random variable. A random variable can be defined
as the numerical quantity associated with the outcome of a random
phenomenon (or an experiment) in such a way that, as the various events
or possible outcomes of the phenomenon occur, the random variable
takes on definite values (Parzen, 1962; Feller, 1968). As examples
the random variable could be the instantaneous point value of the con
centration in a reactor, the instantaneous value of the velocity in a
turbulent flow field, or the displacement of a particle due to turbu
lent diffusion.
3.1.1 The FokkerPlanck equation
Consider a continuous stochastic process X.. As we have seen
1
above, at any time t, X.(t) can be viewed as the numerical value of
a random variable or we can regard X.(t) as the random variable. For
1
a complete probabilistic description of the stochastic process, all
the moments (infinite in number) are required. Alternatively one
should know all the joint probability density functions; i.e. for any
n, p(Xil,t ; X i,t2; .... X. ,t ) should be known. Using the defi
i 1 i 2. 1 n
nition of the conditional probability density this can be written as:
p(X ,t; X. t2; ,t ) = p(X ,t2; ....X ,xn X ,tl )p(X ,tl)
1 2 n 2 n 1 1
(3.1)
Since we are interested in p(X. ,t ), the probability density of the
1 n
n
stochastic process X. at time t we can write Eq. (3.1) as:
1 n
p(X ,tl; X 't2 .... ; Xi ,) =
1 2 n
P(Xi ntlXi tn n; ....; X iti)
n n1 11
P(Xil,t ; X. ,t ); .... X ,t ) (3.2)
1 2 n1
A Markov process can be defined as a stochastic process where
the present state of the process is completely determined by its pre
vious state and is independent of the way in which the previous state
had developed. For a Markov process, the conditional probability
density on the right hand side of Eq. (3.2) can be written as:
P(Xi tli ,tn; Xi 't ; ....; Xit = P(Xin Xi ',tn
n n1 n2 1 n n1
(3.3)
for t1 < t2 < .... t
1 2 n
Eq. (3.3) is the mathematical form of the definition of a
Markov process (Seinfeld and Lapidus, 1974). The conditional prob
ability p(Xi ,t X. ,tn l) is called the one step transition prob
n n1
ability. The probabilistic description of a Markov process does not
require an infinite number of the joint probability density functions
and the problem is reduced to finding the transition probability.
From the definition of a Markov process, it is easy to see that the
transition probabilities obey the ChapmanKolmogorov equation (Bhar
uchaReid, 1960):
p(X,t Xi ,t ) = fP(Xi,tlX!,t')p(X!,t' IX ,to)dX' (3.4)
o o
00
where X! = X.(t') is an intermediate state of the stochastic process
1 1
at t', t > t' > t and where X. and t represent the initial values
S1 0
o
of the displacement of the particle and of the time respectively.
The time t' is fixed in the integration of Equation (3.4). The equa
tion implies that the process can be expressed as going from time t
o
to t in two stages.
Now let us consider an increment of time At. If the present
time in the ChapmanKolmogorov equation is (t+At), if the intermediate
time is t and if the state of the stochastic process (the particle
displacement is (X.AX.), Eq. (3.4) can be written as
1 1
00
p(X ,t+AtIX. ,t ) =j P(X.,t+AtlXnAX.,t)
00
p(XiAX.i,t X ,t )d(AX.)
o
(3.5)
Expanding the integrand on the right hand side of the equation in a
Taylor series in the xdirection about the point p(X.+AX.,t+AtX.,t)
p(Xi,tJXi ,t ) we obtain:
o
p(X ,t+At iX.AX.,t)p(X.AX.,tX. ,t ) =
0
= p(Xi+AXi,t+Atxi.,t)p(X.tIXi ,to) +
0
(_)n
n! ) n
n
J>1
ax"
i
p(XitlXi ,to)]
o
(3.7)
Substituting Eq. (3.7) into Eq. (3.5), we obtain
S
o
t ) = f p(X.+AXi,t+AtlXi,t)p(Xi,tlX,t )d(AX.)
O 1 1
+ j
o n=l
(_)n (A)n
n! i
an
axn
1
[p(X +AX.,it+At Ixi.,t)p.
(Xi,tjXi ,to)]d(AXi)
o
(3.8)
[P(Xi L+!Xi't+AtlXi't)
P(Xillt+Atlxi
= p(X.,t Xi ,t ) / p(X.+AX.,t+AtlX,t)d(AX.) +
o
CO
n
[P(Xi,t Xi ,to)
X [ o
l
(AXi)np(X.+AXi,t+AtjX,t)d(AX.)]
1 1 1 1
(3.9)
From the definition of the probability density, we have
0O
f p(X.+AX.,t+At[X.,t)d(AX.) = 1.
Using Equation (3.10), Equation (3.9) can be rewritten as:
Using Equation (3.10), Equation (3.9) can be rewritten as:
p(Xi,t+AtIXi ,to) = p(Xi,tlX ,t ) +
o o
j(AXi)np (Xi+AX, t+At Ii,t)d(AXi )]
00
n
[p(Xi,tlX. ,to)
i0 o
0
(3.11)
Dividing Equation (3.11) by At and taking the limit as At 0:
Lim
At 0
p(Xi,t+AtjXi ,t ) p(X ,tXi ,t ) =
o 1 o
i (1)n
n!
n=1
an
8 n
X
1
[p(Xi,tli ,t )
0
n
S(AX.)
00
p(Xi.+AX., t+At ., t)d(AX.)
(3.10)
Lim
At + 0
(_l)n
I n!
n=l
(3.12)
The left hand side of the above equation is simply the definition
of the partial derivative of the transition probability with re
spect to time.
3p(Xi,tlxi ,t )
0 
at = Lim
At*O
p(Xi,t+AtIXi ,to)p(X.,tlXi ,t )
o o
th
The n moment of a random variable X. with the probability density
function (X) is defined as (Tucker, 1967)
function p(Xi) is defined as (Tucker, 1967)
E {X} = f Xp(X.)d(Xi)
(3.14)
th
Thus the rate of change of the n moment of the increment AX.
conditioned on (Xt) be defined as
conditioned on (X.,t) can be defined as
1
a. (Xi t) Lim 1 (AX) np(X +AXit+At X.,t)d(AX.)
n 1 At 1 1 1 1 1
1 At0
(3.15)
Substituting Equations (3.13) and (3.15) into equation (3.12), the
"forward FokkerPlanck equation" is obtained.
ap(xi,txi. ,to)
at
n Sai [ac (Xi,t)p(X.,ti x. ,t )]
n! a ni 1 n 0 0
S< X. < t
t > t
(3.13)
(3.16)
The boundary conditions are:
Lim p(X, tI X ,t ) = 6(X.X )
t0 o o
(3.17)
Lim P(X.,tlX. ,t ) = 0
Ixl 1 o
Equation (3.16) is also known as the forward Kolmogorov equation.
The threedimensional FokkerPlanck equation is similar and the de
rivation is identical to the one described above except for Equation
(3.6) which changes to a threedimensional Taylor series and has ad
ditional terms.
3.1.2 Stochastic differential equation
Kirmse (1964) in his Ph.D. thesis describes the displacement
of a particle X. at the end of the time step At as:
im
m
At2
X. = X. + V At + a. m (3.18)
m m1 m1 m1
where X. is the displacement of the particle at the end of the
Im1
st
(ml) step. V. and a. are the velocity and the acceleration re
1 1
th
spectively. The displacement of the particle at the end of the m
step is determined completely from its position at the end of the
previous time step. The single most recent value of the displacement
contains as much information about evaluating X. at t as does the
1 m
m
entire past history of the particle. Therefore the displacement of
a particle in the KirmseFahien model is a Markov process. A sto
chastic differential equation for the model is derived.in this sec
tion considering the displacement of a particle as a Markov process.
A general stochastic differential equation for a Markov
process can be written as (Stratonovich, 1968)
dX. dg
S= f (Xit) + (Xi,t) (3.19)
dt i dt
where X. is a stochastic process with Markovian properties, f.(X.,t)
is the deterministic function and g.(X.,t) is the random function of
the process.
The first two moments of the displacement AX. can be calcu
1
lated as:
a. Li= A 1 f (AX.)p(X+AX.t+Atl.,t)d(AX.)
1. Atmo At f i
10
= E Xi.t = fi(Xi,t) (3.20)
and
a = (dX) Xt (3.21)
21 dt
The deterministic function f. is then the deterministic
1
th
velocity in the i direction in a turbulent flow and the random
function dg./dt is the random velocity contribution in the th
direction. Differentiating Eq. (3.18) with respect to time and
using vector notation, we have
d2X df d2
=+ (3.22)
dt dt
where X is the displacement vector of a particle in space,
d2X/dt2 is its second derivative with respect to time as it moves
in the space and is therefore the Lagrangian acceleration vector.
The derivative df/dt represents the mean Lagrangian acceleration on
2 2
the particle and (d2g/dt2) is the fluctuating Lagrangian acceleration.
Then
df
a (3.23)
d2
S= a' (3.24)
dt2
where a and a' are the mean and the fluctuating components of the
Lagrangian acceleration vector respectively. Equation (3.22) can
now be written as
2
d2X
= a = a + a' (3.25)
dt2
The KirmseFahien model uses empirical probability distri
butions for each of the components of the Lagrangian acceleration
and simulates the path of fluid particles on a digital computer.
The model is discussed in detail in Section 4.1.
The integrated form of Eq. (3.25) is the stochastic differ
ential for a model considering the displacement as a Markov process.
Integrating Eq. (3.25), the stochastic differential equation is ex
pressed in terms of the Lagrangian accelerations and in the Lagrangian
coordinate system as:
dX t t
d= a dt + V' + a' dt (3.26)
dt = o 0 f
Lag t t
o o
where V' is the turbulent velocity vector possessed by the particle
o
at the beginning of the simulation. The above equation in the Euler
ian frame of reference is simply the following:
dX t t
= V + a dt + V' + a' dt (3.27)
dt o J
Eul t t
0 0
where V is the mean velocity vector.
3.1.3 Comparison of the FokkerPlanck equation with the diffusion
equation
The FokkerPlanck equation for a = 0 for n > 3 reduces to
n.
i
S1 a2
S(Xti t) = ( p) + 2 (a p) (3.28)
at 0 3X. 1. 2 2 2.
o l l i
1
th
where X. is the particle displacement in the i direction. The
1
th
first moment in the i direction a from Equations (3.20) and
1.
(3.27) then is
t
= fi(Xit) = V. + f a. dt (3.29)
1
t
o
The first moment in the KirmseFahien model a. is therefore
1
the deterministic velocity which is the sum of the average velocity
in the flow field and another term arising due to the mean Lagrangian
acceleration.
The physical meaning of a2. is easy to see from Equation
1
(3.21), which can be alternatively written as
1 dt 2 (3.30)
a = I(d{ X (3.30)
2 2
where dX. is the expected value of (dX.) conditioned at (X.,t).
1 1 1
38
G.I. Taylor (1921) expresses diffusivity as the time derivative of
the variance of dispersion
1 d 2
i 2 dt i (3.31)
Therefore 1/2 a2. is a point eddy mass diffusivity which is
21
I
a function of the particle position and time. The diffusivity 1/2 a2.
is based upon the Lagrangian statistics of the particles (Lamb and
Seinfeld, 1973) and may therefore be called as the Lagrangian diffu
sivity, KL..
1i
The FokkerPlanck equation can be written in terms of the con
centration instead of the transition probability. The displacement
X. corresponds to the coordinate x, when we go from probability space
1 1
to time averaged concentration space. Then
p(Xi,t) = c(xi,t) (3.32)
t
al (Xit) = V. (xit) + f a (xi,t) dt (3.33)
t
o
Substituting Equations (3.32) and (3.33) into Equation (3.28) and
using 1/2 a2. = KL we have
1 1
39
t
ac (x ,t) 
('t) = [V.(xi.t) + a (x.,t) dt] c (xit)
t
o
2
+2 [KL (x.,t) c(x.,t)] (3.34)
8x 1
The timeaveraged diffusion equation, including the convec
tive flow term, is
ac(xi,t)
at 3 x. [Vi(x.,t) c(xi't)]
+ [K (x ,t) c(xi,t)] (3.35)
ax. i i' Dx. 1
1 1
The Eulerian diffusivity K. can be related to the Lagrangian diffu
1
sivity KL by equating Equation (3.34) and (3.35) and solving for
1
t
c a dt +a (L c)
1 1
t
K. o) (3.36)
1 (ac/Ix)
where t represents the conditions at the source. The two diffu
sivities are equal to each other if homogeneous conditions exist and
if the Lagrangian accelerations are equal to zero.
The Lagrangian diffusivity in general should be expected to
be a function of the source and of the initial time as well as of
the time of diffusion. A question then arises as to whether the
Eulerian diffusivity K. should also be a function of the source.
1
Molecular diffusivity is a property of the medium and if an eddy dif
fusivity is to be defined in an analogous manner, it should be a
property of the flow field only rather than a quantity associated
with the source and the time of diffusion.
The rate of change of the moments of the third and higher
order a 's are neglected to obtain correspondence with the diffusion
ni
equation. If the diffusion process has Markovian properties, then
in general the moments of higher orders should be nonzero and should
be retained in the FokkerPlanck equation. The Boussinesq hypothesis
of the flux being proportional to the concentration gradient in that
case is invalid. Several researchers in the field of turbulence have
questioned the basic concept of a turbulent diffusivity, and this
criticism does seem to have some validity. However the diffusivities
are useful empirical coefficients in the flux laws and possess many
advantages in routine dispersion calculations.
3.2 The Growth Law
3.2.1 The Lagrangian similarity theory for large times
The concept of Lagrangian similarity was first introduced by
Batchelor in 1957. If homogeneity exists, a particle moving in the
downstream direction finds itself surrounded by the same turbulent
structure. For free shear, steady state flows, it is assumed that
the structure of the flow is the same at different downstream dis
tances. Mathematically, this assumption is expressed as:
u(t) u
(t) = F() (3.37)
W(t)
where u(t) is the instantaneous velocity of the particle, u is its
mean velocity, w(t) is a velocity scale of the particle motion at
time t, and the independent variable n is proportional to the ratio
of the time dt to T the time scale of the particle motion.
dt
dr a d (3.38)
T(t)
The function F is stationary and random. The statistical
mean velocity of all the particles is then obtained by averaging Eq.
(3.37) over all the particles. Thus:
dX
dt = u + F(n) w(t) (3.39)
where X is the displacement of the particle. The function F(n) is
obtained by assuming that the length and the velocity scales of the
particle motion are given respectively by:
_ P
(3.40)
2,(t) t [X(t)] 
and
W(t) a [X(t)] 2 (3.41)
The time scale of the motion can be obtained by taking the
ratio of the length to the velocity scale.
Pl+P2
T(t) [X(t)] 1 (3.42)
Batchelor (1959) later suggested that one should be able to obtain
the same results for a ground level steady state source in the atmos
pheric ground layer. He developed this idea in an unpublished work
which is discussed by Ellison (1959), Cermak (1963) and Csanady (1973).
It was assumed that the Eulerian properties of the flow structure de
pend only upon the friction velocity in a neutrally stratified layer
and that the roughness length affects the mean velocity of the flow
only. The similarity hypothesis for elevated sources is expressed
by Cermak (1963) as follows:
For a marked particle which is at z = H when
t = 0, the statistical properties of particle motion
at time t depend only upon u* and ttv where t is of
order H/u* or larger, where t is a virtual time ori
gin with magnitude of order H/u*. (page 51)
The mean longitudinal and the vertical positions of the par
ticles can be expressed as
dZ
dt= bu (3.43)
dX t ) (3.
dt z
o
for t>> H
U*
The concentration is assumed to be of the form
m1
c a x (3.
and the plume width at the ground is assumed to be of the form
Y m2 (3.
Y a x (3.2
44)
45)
6)
Based on the two above mentioned relationships, the exponent for
the plume width at the ground is found to be
x
z
m2 log 
z Z
0 0
(3.47)
However, Batchelor's hypothesis (1959) and Cermak's (1963)
developments are for a time of releases larger than H/u*. For a
stack height of 200 m and a friction velocity of 0.4 mps, t should
be larger than 500 seconds. At a velocity of 6 mps, one must be at
and
least 3 km away from the source. Cermak (1963) compared the theory
with some atmospheric diffusion experiments and with some wind tun
nel data. He concluded that the ratio of source height to roughness
length H/zo is of major importance in atmospheric dispersion calcu
lations.
3.2.2 The growth law for a Markov process
The Lagrangian similarity theory of Batchelor is modified
for small times based on the correlation coefficient results.
For a Markov process, the correlation coefficient is of the
exponential decay type (Csanady, 1973).
T
RL(T) = e (3.48)
where T is the integral Lagrangian time scale. If Z represents the
displacement of a particle in zdirection, we can write (Taylor, 1921)
2 2 t tt T
22 2 d 2
 z = 2 w'2 e dT dt' (3.49)
dt dt2 0
0 0
2
where a is the variance of the dispersion in the zdirection.
z
For neutral conditions (Panofsky and McCormick, 1960; Cramer
and Record, 1969)
w2 u2 (3.50)
w2 2
w' = blu
1 *
where bl is a universal constant.
Let x, y and z be the mean positions of the particles re
leased from the source. Since the rate of change of z is very small
compared to the rate of change of x, and since axial diffusion can
dx dx
be neglected, one can write = as the mean velocity of the
dt dt
fluid at z = z (Cermak, 1963). For a logarithmic velocity profile
dx u* z
dt K Z
d~ir11 ^"o
(3.51)
Substituting Equations (3.50) and (3.51) into Eq. (3.49), we have
K* z
( 0 
d22
dZ 2
2 = 2blu
dx2
X Xt X"
Sf e L dx" dx' (3.52)
0 0
Integrating Eq. (3.52) twice with respect to x, we obtain
the following expression:
2
2 2b K
z = [in +i 2
[lnZ +1)]2
z
o
a =bA 1+
z L AL l
Sx 
2 x AL
A 2 A 1 + e L
L AL l
L
2 3
x x x
A1 +
L2 2!nZ 3!AL3
2 2b K
b 2
In( F + 1)1
where
(3.53)
1
2
(3.54)
(3.55)
Eq. (3.53) can be written as:
A x [ 2 3 1
L x x x 2
z = b 1 + + ... 2
z /TA 3A 12AL A3
L L
(3.56)
Since z changes very little with x, it may be replaced by the stack
height H. Eq. (3.56) can be nondimensionalized by H to obtain an
*
expression for the dispersion coefficient a = az/H:
z
S2 3 
b x x 2
S 1 ++(3.57)
S 12where
where
/ 2b K
b 1
ln(+ 1
*
z
z 
o H
(3.58)
(3.59)
Expanding Eq. (3.57) we obtain:
*
b x
S= x* 1 
z r 6AL
*2 *3
x x
x x
36AL 270AL
The infinite series above can be approximated by:
(3.60)
47
0 = x (3.61)
z 1
2 1 + x
6AL
or
oz =z (* (3.62)
1 1 + b x
The growth law for a Markov process should therefore be given
by Eq. (3.62). It is important to see if this behavior is obtained
for the empirical correlation coefficient curves as well.
Since the Lagrangian correlation coefficients are assumed to
be similar to the Eulerian correlation coefficients (Sec. 4.2), it is
possible to obtain the shapes of z vs. t curves by integrating Eq.
z
(3.49) twice.
St t'
O a RE(T) dT dt' I (3.63)
S0 0
The experimental data of Prairie Grass was used for the cor
relation coefficient to obtain the right hand side of Eq. (3.63)
which we will call It. The resulting curves are fitted to the form
t
I+ b (3.64)
t 1 + b2t
This curve fit was carried out for all the stability condi
tions. The results for one typical neutral (Run #24), one unstable
(Run #61) and one stable (Run #22) conditions are shown in Figures
(3.1) to (3.3). The shape of the growth law obtained from the cor
relation curves is in excellent agreement with Eq. (3.64) which is
obtained by assuming that the diffusion process is a Markov process
and consequently the correlation curves are of an exponential decay
type.
The same form of the growth law is expected to be valid for
the dispersion in the ydirection as well. Thus
*a bx (3.65)
y H 1 + b x
Y2
Although the effect of stability conditions on the form of
the growth law is assumed negligible, the parameters are functions
of stability conditions, and they should be correlated with the proper
stability parameters.
 / Experimental
t, sec
Figure 31. Comparison of It curve calculated from the correlation coefficient data
with the Markov process curve fit for typical neutral conditions.
24
H 16
0
Experimental
8
Markov process curve fit
0 II I I.
0 10 20 30 40 50 60 70
t, sec
Figure 32. Comparison of I curve calculated from the correlation coefficient data
with the Markov process curve fit for typical unstable conditions.
1.6
1.2
0.8
0.4
0
t, sec
Figure 33.
Comparison of It curve calculated from the correlation coefficient data
with the Markov process curve fit for typical stable conditions.
CHAPTER 4
THE STOCHASTIC MODEL
4.1 Description of the Model
Diffusion from a steady state, continuous point source is
modeled in this work. For the application of the model used in this
work, the origin is at the ground level directly below the source,
the xaxis extends horizontally in the direction of the mean wind,
the yaxis is in the horizontal plane perpendicular to the xaxis
(cross wind), and the zaxis is in the vertical direction. The con
tinuous emission of pollutants from a stack is assumed to be composed
of a finite large number of mass pockets (of the order of 104). Each
of these mass pockets is referred to as a particle. Each particle
is released from the source with a statistically random velocity.
The motion of the particle in accordance with a stochastic equation
is simulated by a Monte Carlo method until it exits the region of
interest.
Statistical analysis of a velocity field in a steady state
turbulent flow indicates that the velocity has an approximate mar
ginal Gaussian distribution with constant variance, and that the
time derivatives of the velocities are also continuous and also have
an approximate marginal normal distribution (Simmons and Salter, 1938;
Townsend, 1947). Gaussian distributions with measured or correlated
variances are used to simulate the velocities of the fluid particles.
The variances of the accelerations are calculated from information
regarding the velocity correlations and the relations between the
Lagrangian and Eulerian correlations.
For a digital simulation, it is necessary to discretize the
probability distribution curves (See Appendix A). The particles are
moved in finite time steps, but the simulated instantaneous velocity
field is continuous. In the previous chapter, an equation for the
instantaneous acceleration was derived (Eq. 3.27) from the general
stochastic equation for a Markov process. The displacement AXim of
th th
the particle in the i direction during the m time step At is ob
mt
trained by integrating Eq. (3.25) twice in the i direction over At .
m
Then
AX, = V At + a. dt + a! dt, (4.1)
m i m i I
At At
m m
where V. is the total velocity of the particle at the beginning
m1
of the time step. If the Lagrangian. accelerations are assumed con
stant during the time step At ,
At 2 At 2
m m
AX. = V. At + a + a! (4.2)
1 1 m i 2 1i 22)
m m1 m1 ml
The total velocity V. is the sum of the mean velocity V. and the
i i
m1 m1
fluctuating velocity. The model is capable of using any given form
of the mean velocity V. (x,y,z). The mean Lagrangian accelera
m1
tions are obtained from the equation of motion by estimating dif
ferent terms from the measured meteorological parameters. The fluc
tuating accelerations are obtained randomly from Gaussian distribu
tions (See Appendix A). The time steps are random and without any
memory, and thus can be obtained from an exponentially distributed
population (Kirmse, 1964). The mean time step can be looked upon as
a time scale of turbulence.
Since it is assumed that no interaction takes place between
individual particles, they are released from the source in succession.
The position of the particle at a fixed interval of time t is stored
s
in the computer to get the concentration values. The concentration
is calculated from the particle sampling according to the following
equation.
Q Ncell(ts)
c= Q cell (4.3)
N V
p cell
where
1
Q = source strength, MT
N = number of particles released
N = number of particles positions
(samplings) in a cell
t = sampling time, T
V = volume of a cell, M
cell
It can be shown that
Q is the source strength represented by each particle.
N
p
Qts
 is the mass of the pollutant represented by each sampling.
N
p
Qt Ncell
s cel is the mass of the pollutant given by all the samplings
N
p in a cell.
Total reflection is assumed at the ground; i.e., when the
particles hit the ground, they are reflected in a perfectly elastic
manner. The molecular motion may be neglected in the atmosphere be
cause of very large time and length scales of turbulence (Monin and
Yaglom, 1971; Tennekes and Lumley, 1972).
It is assumed that as soon as a particle is released it
adopts the flow properties of the environment. If a temperature or
a velocity difference exists between the fluid coming out of the
stack and the surrounding fluid, their flow properties will be differ
ent, and thus the source height should be taken as the effective
stack height and not as the physical height.
4.2 Estimation of Simulation Parameters
A particle is followed as it moves downwind and therefore one
needs the Lagrangian turbulent properties. However because of the
fixed coordinate system used, the mean velocity profile must be known.
For the simulation of Eq. (4.2), the following parameters are
needed: the standard deviations of each of the components of the
fluctuating velocity (turbulence density), the mean Lagrangian ac
celeration profile, the standard deviations of the fluctuating com
ponents of the Lagrangian accelerations, and the mean values of the
time steps.
Since only the Eulerian properties are usually measured and
reported, Lagrangian properties are obtained through the use of the
factor
B E A/AL (4.4)
where AE and AL are the Eulerian and Lagrangian macroscopic length
scales respectively. Kirmse (1964) and Gielow (1972) used a value
of.B = 0.7. Hay and Pasquill (1959) propose a constant ratio, g of
the Lagrangian and Eulerian time scales, i.e.
RL(St) = RE(t) (4.5)
where RL and RE are the Lagrangian and Eulerian velocity correla
tion coefficients respectively. Their estimate of B varied from 1.1
to 8.5 Angell, et al. (1971) reported B from 2.5 to 4.5 However
they estimated B from the ratios of the frequencies at which the
maximas are observed for the Lagrangian and the Eulerian velocity
spectra. These are point values of 8 and do not represent the values
that might be satisfactory over the complete range of frequency
(for spectra) or time (for correlations). Tslitzer and Dumbauld
(1963) estimated B = 5 for all their runs. Since their value is
intermediate between those of Hay and Pasquill, this value is used
throughout this work.
The factor B is calculated from B as follows:
Let TL and TE be the Lagrangian and Eulerian macroscopic
time scales respectively. Then,
TL =J R (Bt) d(Bt) (4.6)
0
and
TE = RE (t) dt (4.7)
0
Expressing TL in terms of TE
TL =fR (st) d(St)
0
j RL (st) dt (4.8)
0
Using Eq. (4.5)
TL = RE (t) dt = E Tg (4.9)
0
= TL/TE
B = AE/AL
EL \\
(4.10)
(4.11)
Using well known relations between the space and time macroscales
(Hinze, 1959)
AL v TL
AE = u TE
and therefore
A uT 
E UTE u
B=
A T 6o
V V. v.
i L
th
where a is the intensity of turbulence in the i direction.
V.
1
(4.12)
(4.13)
(4.14)
By taking the Fourier transform of Eq. (4.5), the Lagrangian
and Eulerian energy spectra can be correlated as follows (Hay and
Pasquill, 1959):
FL(n) = 6 FE(n)
I
(4.15)
where FL and FE are the Lagrangian and Eulerian energy spectra
respectively, and n is the frequency. Eq. (4.15) when integrated
over the entire frequency range gives the equality between the
Lagrangian and Eulerian turbulence energies. But Bullin and
Dukler (1974) used 0 = 4 in the following equation which is an
approximation of Eq. (4.15).
FL(n) = FE(n) (4.16)
A value of 8 = 4 would imply that the Lagrangian turbulence
energy is four times as large as the Eulerian turbulence energy and
this would lead to serious error in using their model.
The other parameters are calculated as follows:
4.2.1 Neutral conditions: Complete homogeneity is assumed in both
the x and the y directions.
4.2.1.1 The mean velocity profile
The velocity profile for neutral conditions is of the loga
rithmic form (Monin and Obukhov, 1954). Many derivations are avail
able in the literature, and the result is independent of the method
of derivation (Lumley and Panofsky, 1964). The derivation presented
here is based on dimensional analysis.
The roughness of an atmospheric ground layer is described
by a roughness length z (also known as "dynamic roughness"). The
roughness length is not equal to the actual height of the vegetation
60
or other physical structures on the ground, but is considerably
smaller in magnitude. It varies from the fraction of a mm over
snow or sand to several meters over cities (Nawrocki and Papa, 1963).
Since the mixing length in the atmosphere is proportional to
(z+zo), the characteristic length scale is taken as (z+z ). The
characteristic velocity in the atmosphere is the friction velocity.
u = (4.17)
where T is the shear stress at the ground and p is the density of
o
the atmosphere. The frictional velocity represents the drag on the
surface. By dimensional analysis, the velocity gradient ( /dz) is
therefore proportional to u /(z+z ). Taking the constant of propor
0.
1
tionality as , we have
I u
du 1
(4.18)
dz K (z+z )18)
o
Integrating from z = 0 to z = z, we obtain
u z+zo
u = in o (4.19)
K z0
The dimensionless constant K is the von KarmAn constant, and is taken
as 0.4.
The roughness length z is a property of the layer and is
L
61
independent of meteorological conditions. The friction velocity
u* represents the magnitude of the wind velocity and is invariant
with height.
4.2.1.2 The variance of the velocities
These are assumed to be equal to the Eulerian mean square
fluctuations.
o2 = ( ) = (a= 2 (4.20)
v. v. v. E
4.2.1.3 The variance of the accelerations
The Lagrangian acceleration is obtained from the Eulerian
acceleration by the LangrangianEulerian transformation using the
ratio of the length scales B (Kirmse, 1964).
 2
a =B a ((u)2 (4.21)
a. v. x
ai 1
The standard deviations of the velocity can be obtained either
from the experimental data or from correlated data. The correlation
of the spatial derivative of the velocity can be obtained from the
energy dissipation e if one assumes that small scale structure of
the turbulence is isotropic. This assumption is very reasonable be
cause isotropy is being assumed only at the dissipation scale of tur
bulence while the overall behavior is still anisotropic. Then
(Hinze, 1959),
2
E = 15V' ( ) (4.22)
The energy dissipation has been found to be a function of
the wind velocity and height (Record and Cramer, 1966). Record and
Cramer added their Round Hill data (Cramer et al., 1961) to the
composite graph of Ivanov (1962). The dissipation rate and its
height dependence were found to vary significantly with thermal strat
ification.
In a later report, Cramer and Record (1969) were able to suc
cessfully correlate all their energy spectra measurements by plotting
F(n)u
all of their data on a graph of dimensionless spectra 2 vs.
u* z
nz
dimensionless frequency This relationship gives the height
and the velocity dependence of the energy spectra, and consequently
the energy dissipation rate. The latter is related to the energy
spectra by the following relation:
0
e = 2v/ k2E(k) dk (4.23)
o
where k is the wave number.
2k n
k = (4.24)
_9
L
E (k) =  F(n)
2Tr
(4.25)
The wave number k can be nondimensionalized in a similar manner to
the frequency n if we let
* nz
n 
u
(4.26)
Multiplying Eq. (4.24) by z and substituting Eq. (4.26) into
it, we have
kz nz *
2  n
2f d
(4.27)
2
Similarly, dividing Eq. (4.25) by (u*z) and simplifying, gives
27TE u F *
=~ =~ F
2 2
u*z u*z
(4.28)
*
Substituting for E and k in terms of F and n in Eq. (4.23), we
obtain
(4.29)
0o + 2 2
E = 2v J 2n*z ) u z d ( 2
0
and
64
872r 2 2 ,2 ,
E 82 n F dn (4.30)
z
0
2 2
87r2 u* *
2 (4.31)
z
where
2 2
n F dn = (4.32)
0 87r v u
0
The dimensional energy dissipation rate c is taken as a
*
constant and its value obtained from F vs n curves (Cramer and
Record, 1969). Combining Equations (4.21), (4.22) and (4.32), we
obtain an expression for the standard deviation of the Lagrangian
accelerations.
S= 2ru* B /2E* (4.33)
ai z vi
4.2.1.4 The mean time steps
Kirmse (1964) and Gielow (1972) have shown that the mean
value of the time step is proportional to the ratio of the standard
deviations of the velocity and of the acceleration, or
a
v.
X = 0.5 (4.34)
t. a
1 a.
1
1
where the factor 0.5 is determined theoretically.
4.2.1.5 The mean Lagrangian accelerations
The mean Lagrangian acceleration in the xdirection a is
x
the derivative of the turbulent shear stress (Kirmse, 1964; Gielow,
1972). The variation of the shear stress with height is extremely
small in the atmosphere (Monin and Obukhov, 1954; Lumley and Panofsky,
1964) and a is therefore taken as zero. Because of no net buoyancy,
x
the mean Langrangian acceleration in the zdirection is also taken
as zero.
4.2.2 Nonneutral conditions
4.2.2.1 The mean velocity profile
In the case of nonneutral conditions, one must use an addi
tional length scale called the MoninObukhov length L. The expression
for the shear then changes to
= Y ^(4.35)
dz K(Z+Z) (4.35)
where I is a universal function of Z/L. Since L = for neutral con
z+z
o
editions, = 0 and therefore
L
p (0) = 1
(4.36)
66
The dimensionless shear i can be expanded in a Taylor ser
ies about zero and for small values of /L (< 0.1), one can
write
z+z \ z+z
L = 1 + B  (4.37)
Substituting Eq. (4.37) into Eq. (4.35) and integrating from z = 0
to z = z, we obtain
u = [ in ( z + B (4.38)
K z L
The MoninObukhov length is defined as follows in terms of
the turbulent heat flux:
3
U*
L = ) (4.39)
K(g/TR) (q/CpP)
Since the heat flux is difficult to obtain in practice, a convenient
way to avoid this difficulty is by expressing the momentum and heat
fluxes in terms of the velocity and potential temperature gradients.
SKm du (4.40)
pC dz
PC
67
where K and Kh are the momentum and heat turbulent diffusivities
m
respectively. Assuming that the variation in T with height is very
small and can be neglected, we have
T 
o 2 du
= u = K d
p m dz
(4.42)
Substituting Equations (4.41) and (4.42), Eq. (4.39) becomes
du
uL (Km dz
L=
R
K u, (du dz)
T R idz
(4.43)
(4.44)
A gradient
similarity length LG is defined as
(du
LKh dz
L LL
G Km K dz
m T R dz
Then the velocity profile equation can be written as
u z+z
u In + G
K z\o) G LG
SG
for zIL < 0.1
L G
G
(4.45)
(4.46)
where BG is a universal constant. The constant G was taken as
4.5 for unstable conditions and 7.0 for stable conditions (Lumley
and Panofsky, 1964). For IZ/LI > 0.1, Equation (4.37) is not valid,
but a power law equation (Csanady, 1973) can be used:
 u (zl
u 1 (4.47)
where zl is an arbitrary height. The exponent a depends upon the
roughness length zo, the reference height zl and the stability.
Since the loglinear profile given by Eq. (4.46) can be assumed to
hold up to the height IZ/LG I the reference height
z1 = (4.48)
was chosen as this value. Therefore,the reference velocity is
LG LG
u G+ z
u = L* In( 10 z) +S G 1(4.49)
Sz L
o G
= (zo LG) (4.50)
The gradient similarity length LG describes the stability while zo
describes the physical characteristics. The only dimensionless param
eters that can be formed from these parameters is Zo/LG. Therefore
1
S= a L 0 (4.51)
The above equation is seen to be obeyed for all the Project
Prairie Grass runs (Barad, 1959). The exponent a is plotted against
the stability parameter z/LG in Fig. (4.1) for all the runs with
stable conditions. The function a can be correlated for stable con
ditions by the expression
a= a In ) + b (4.52)
Zo
where a is the slope ofa vs. ln( /LG) line and b is the intercept
ZO
at /LG = 1. For unstable conditions
a = a constant
= 0.124 (4.53)
4.2.2.2 The temperature profile
Potential temperature is defined as the temperature a parcel
of air would attain if it was brought to a standard reference pressure
under adiabatic conditions. The potential temperature can be related
to the temperature in the following manner:
Assuming the air to behave as an ideal gas, the first law of
thermodynamics for adiabatic conditions can be written as (Hess, 1954)
RT
C dT dP = 0 (4.54)
p p
0.7
0.6
0.5
0.4
0.3
0.2 
0 0
Q.1
0 2
105 22
10 2
O0
OcD
0
2
zo
ZO/L5
5 102
Figure 41. The velocityprofile exponent a vs. the stability parameter z** for stable
conditions.
I r 1
Integrating Eq. (4.54) from the reference pressure PR to P, we have
PR
S= T () Cp (4.55)
where 8 is the potential temperature. For the Prairie Grass runs
(Barad, 1959), the pressure profiles were observed to be linear. In
this work, the reference pressure is taken as the pressure at the
ground. Thus the pressure profile is
P = PR + PR g z (4.56)
where pR is the density at the ground, and gz is the acceleration
due to gravity in the zdirection, which is given by
gz = g = 9.8066 m/sec2 (4.57)
Substituting Eq. (4.56) into Eq. (4.55) and expanding into a bi
nomial series, we have
PR gZ z R
9 = T(1 ....) (4.58)
PR C
Using the ideal gas law and neglecting the higher order terms, Eq.
(4.58) reduces to
0 = T +  z (4.59)
C
p
The above equation is a good approximation for Eq. (4.55)
for small values of z.
The potential temperature 6 was plotted against the mean
velocity u (Figs. 4.2 and 4.3) for a large number of Prairie Grass
stable and unstable runs. The slopes of these curves were observed
to be constant. This implies that the gradient similarity length
LG can be taken as a constant in the surface layer (Eq. 4.45). The
potential temperaturemean velocity curves can be expressed as
0 = R + u) (4.60)
Combining Equations (4.46) and (4.60) we have
6 = 6 + ln z + B (4.61)
R du K z GLG
A gradient temperature scale TG analogous to the friction velocity
is defined as
dO
T dz u, (4.62)
G di K
dz
The loglinear temperature profile equation then is
( z+z
e e= R6+T rin 0oo + z (4.63)
R G Lz G LG
10
0 2 4 6 8
u, mps
Figure 42.
The potential temperature vs. the.mean velocity
for some Prairie Grass runs with stable conditions.
74
36
#45
#61
32
#5
30
28 
#16
#152 0
26
24
22 I I I
0 2 4 6 8
u, mps
Figure 43. The potential temperature vs. the mean velocity for
some Prairie Grass runs with unstable conditions.
The gradient temperature scale TG is related to the MoninObukhov
temperature T* according to the following equation.
K
T = T
G Kh
(4.64)
Similarly the power law temperature profile is
T i
S= eR + I
R (0.1)t L G
(4.65)
4.2.2.3 The variance of the velocities
4.2.2.3 The variance of the velocities
As in the case of the neutral conditions, the velocity vari
ances are assumed to be equal to the Eulerian values which may be
obtained from the experimental or correlated data.
4.2.2.4 The variance of the accelerations
The variation of the energy dissipation with stability is not
available. But the functional relationship of with z is assumed
to be the same as for the neutral conditions. This leads to the
0 being inversely proportional to the height. The correlation co
efficients were measured during the Prairie Grass runs (Haugen, 1959),
and therefore the microscales of turbulence could be estimated. From
Hinze (1959), the time microscale tE. is related to the correlation
1
76
of the velocity derivative according to the following equations:
(4.66)
1 1
2 2
tE
1
f(xi) = RE.
1
(4.67)
where f(x.) is the longitudinal space correlation coefficient. Also
1
ax
2 i_2
vi x.
1
(4.68)
X. =0O
1
According to Taylor's hypothesis
t u
3t ax
(4.69)
Therefore
(u' x)2
2
a
V.
= 2
2 2
u tE.
1
Using Equations (4.14) and (4.21), we obtain
r2 a
v.
ai BtE
1 E
(4.70)
(4.71)
4.2.2.5 The mean time steps
The mean time steps are obtained from the turbulent intensity
and the variance of the acceleration in the same manner as for the
case of neutral conditions (Eq. 4.34).
4.2.2.6 The mean .Lagrangian accelerations
The mean accelerations in the x and y directions are taken
as zero. The mean Langrangian acceleration in the zdirection is
obtained from the equation of motion in the zdirection.
Dw 1 3P
a t = t + g (4.72)
z Dt p z z
Let P (z) describe the pressure at steady state neutral conditions.
Then
P = PR + P g z (4.73)
If 3 represents the deviation of the density from po, then
1 1 1 (474)
P + P Po (4.74)
1_ (4.75)
PO po2
0o
78
The deviation p can be calculated by expanding p into a twodimen
sional Taylor series about the point po as follows:
T+ P+ (4.76)
P Pp0 T Po
where
T = T To, (4.77)
P = P P (4.78)
and T and P are the mean temperature and pressure profiles re
o o
spectively for steady state neutral conditions. Assuming the ideal
gas law and neglecting the higher order terms,
P = P + P (4.79)
o o
Substituting Eq. (4.79) into Eq. (4.75), we obtain
1 1 1 T P
po + p (4.80)
P p p2 OT P
1 1 T (4.81)
p p T
0
Substituting Eq. (4.81) into Eq. (4.72), we have
 1 dP T
a = 1 + + gz
z p dz T P
a0
(4.82)
Since the pressure and temperature profile are known one can there
fore obtain an expression for a in terms of z. The temperature pro
z
file for the neutral conditions is given by
T = 68 + r z
o R
(4.83)
where r is the adiabatic lapse rate and is equal to (g/Cp). Since
the reference pressure is taken as the pressure at the ground level
(4.84)
OR = TR
Substituting Equations (4.59), (4.63), (4.65) and (4.83) into Eq.
(4.77), we have the following expression for T:
( z+z \
S= TG i zo G LG'^
LG
10
T T
G
S 0 1)
(0.1)o
a
z
LG
LG
z> 10
10
(4.85)
(4.86)
1
Combining Equations (4.85) and (4.86) with Eq. (4.83), we have
z+z \
G n G +
T o G
(4.87)
P
O
C
p
+ + **+
+R+ N z
R R 1
for z < L
10
LZG _
(4.88)
(4.89)
( I + + 1 **)
RoD +~ R1
LG
for z > LG
10
where
* z
z =
H
z
o
H
** z
z L
LG
T
o
and
T
T
o
(4.90)
(4.91)
(4.92)
8R
TG
(4.93)
(4.94)
FLG
N 
"1
Similarly, the expression for the pressure deviation term becomes
P + dP
R R dz)
dP
z P( z
R dz
(4.95)
PR+ (dz)
J dP
P dz dz
S o
A linear approximation for the above equation is
P
P
SdP
z dP o 0]
PR dz dz
R
(4.96)
(4.97)
(4.98)
N N 
** 2 1 R
NC
1+ R z
R
d 
(4.99)
PPo
p op
0 0
where
82
Substituting Equations (4.57), (4.98) into Equation (4.82), we
obtain
a 1z R 2 1 R (4.100)
Cp 1 L T + R z*
where /T is given by Equations (4.87) and (4.88). All the terms
in Eq. (4.100) are plotted against the height for a set of typical
conditions in Fig. (44). It can be seen that the variation in a
close to the ground is considerable.
close to the ground is considerable.
0.04
0.02
C 1 C **
N N N z
R 2 2 iR
p S1+N **
0 1 R
0.02
R N2
1 " 2
Cp N
0.04
0 250 500 750 1000
z, m
Figure 4.4. The variation of different terms in the vertical mean Lagrangian acceleration
equation with height.
CHAPTER 5
RESULTS
5.1 Computational Parameters
The Monte Carlo simulation of the stochastic model used in
this work can be carried out for any arbitrarily chosen number of
particles, sampling time and cell size. A sensitivity analysis was
carried out studying the effect of these parameters, and optimal
values then were used for all the runs. In addition, one is free to
divide the distribution functions for the velocities, accelerations
and time steps intoanynumber of parts (see Appendix A). Gielow
(1972) tried different numbers of partitions in his study and even
tually used 16 partitions. Since the requirements of precision in
the atmosphere are less than in the laboratory scale simulations,
sixteen partitions were also used in this work for all the distribu
tion functions.
3 5
The number of particles was varied from 10 to 10 As one
may expect, computer time requirements are roughly proportional to
the number of particles released. The scatter observed in the results
decreases with an increase in the number of particles. However the
marginal utility decreases rapidly with the increase in the number of
particles. From these considerations, 10,000 particles were used in
all the runs.
The sampling time has a similar effect on the results as
does the number of particles. An improvement in the results is ob
served if the particles are sampled more often; i.e. if the sampling
time is decreased. From a similar analysis, a sampling time of 0.25
seconds was selected.
The available computer storage is an obvious limitation on
the total number of cells. If the downwind distance that the parti
cles have to travel increases, then the cell size must be increased.
Also since the spread increases approximately linearly with the down
wind distance, the cell size was varied with the downwind distance;
and instead of a complete network of cells, twodimensional cell
structures were used at different downwind positions. Since perfect
symmetry is assumed in the ydirection, all the particle positions
are folded in the ydirection. The axis corresponds to the center
of the first cell, and is therefore one half the size of the other
cells. The cell size has no effect on the results except when an ex
tremely large cell size is used close to the source.
5.2 Validation of the Model
The model was validated by comparing the results of a simu
lation with those obtained in the Prairie Grass diffusion experiments
(Barad, 1958, 1959; Haugen, 1959). The correlated turbulence data
of Cramer and Record (1969) were used to study the effects of differ
ent physical and meteorological parameters including the effects of
neglecting the axial diffusion. The results are correlated in a use
ful form, and graphs for correlated dispersion coefficients and eddy
diffusivities are plotted.
5.2.1 Neutral conditions
Prairie Grass diffusion experiments (Barad, 1958, 1959; Hau
gen, 1959) were conducted at O'Neill, Nebraska. Sulfur dioxide was
released at a height of 0.46 m from the ground and concentrations
were measured in the horizontal direction along five concentric arcs
and in the vertical direction on six towers. The five arcs were semi
circular and of radii 50, 100, 200, 400 and 800 meters. All the six
towers were along the 100 m arc. For the same conditions as those of
Run #45 (which was approximately neutral), the concentration distri
butions in the ydirection at the 50 m, 100 m and 200 m arcs, and in
the zdirection for two towers (y = 5 m and y = 18 m) are compared
in Figures (51) to (55), respectively, with the experimental re
sults and with the simulation results of Bullin and Dukler (1974).
The concentration profiles, obtained from the PasquillGifford equa
tion (Turner, 1967) are also plotted on the same graph using the
standard deviations for the horizontal and vertical dispersion from
the PasquillGifford charts (class D). In all the cases, complete
homogeneity in both the x and the y directions is assumed. It can
be seen that the results are significantly better than using either
of the two above mentioned models.
Run number 24 was seen to be the closest to perfectly neutral
conditions. The results at different downwind distances in the y and
zdirections are plotted in Figs. (56) to (59). The measured values
of the velocity profile and the turbulent intensities were used in
Near neutral conditions
(LG = 125 m)
x = 50 m
z = 1.5 m
I 
0 Experimental
7
O This work
 Bullin & Dukler (1974) 
PasquillGifford
(class D)
0 0
 I O
\
S .0
0/
^Io
\
I 0@
I O
Figure 51.
y, m
Concentration distribution in the ydirection for Project Prairie Grass
Run #45 at x = 50.
500
400
300 I
200 r
to
1 i .
o\
^
1
