Title Page
 Table of Contents
 Previous studies
 The stochastic model
 Conclusions and recommendation...

Group Title: Statistical modeling and simulation of atmospheric turbulence,
Title: Statistical modeling and simulation of atmospheric turbulence
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00084180/00001
 Material Information
Title: Statistical modeling and simulation of atmospheric turbulence
Physical Description: xiii, 222 leaves. : illus. ; 28 cm.
Language: English
Creator: Pahwa, Suresh Badrinath, 1947-
Publication Date: 1975
Subject: Atmospheric turbulence   ( lcsh )
Chemical Engineering thesis Ph. D
Dissertations, Academic -- Chemical Engineering -- UF
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
Thesis: Thesis -- University of Florida.
Bibliography: Bibliography: leaves 206-221.
Statement of Responsibility: by Suresh B. Pahwa.
General Note: Typescript.
General Note: Vita.
 Record Information
Bibliographic ID: UF00084180
Volume ID: VID00001
Source Institution: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: aleph - 000582317
oclc - 14101024
notis - ADB0691

Table of Contents
    Title Page
        Page i
        Page ii
    Table of Contents
        Page iii
        Page iv
        Page v
        Page vi
        Page vii
        Page viii
        Page ix
        Page x
        Page xi
        Page xii
        Page xiii
        Page 1
        Page 2
        Page 3
    Previous studies
        Page 4
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
    The stochastic model
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
        Page 131
        Page 132
        Page 133
        Page 134
        Page 135
        Page 136
        Page 137
        Page 138
        Page 139
        Page 140
        Page 141
        Page 142
        Page 143
        Page 144
        Page 145
        Page 146
        Page 147
        Page 148
        Page 149
        Page 150
        Page 151
        Page 152
        Page 153
        Page 154
        Page 155
        Page 156
        Page 157
        Page 158
    Conclusions and recommendations
        Page 159
        Page 160
        Page 161
        Page 162
        Page 163
        Page 164
        Page 165
        Page 166
        Page 167
        Page 168
        Page 169
        Page 170
        Page 171
        Page 172
        Page 173
        Page 174
        Page 175
        Page 176
        Page 177
        Page 178
        Page 179
        Page 180
        Page 181
        Page 182
        Page 183
        Page 184
        Page 185
        Page 186
        Page 187
        Page 188
        Page 189
        Page 190
        Page 191
        Page 192
        Page 193
        Page 194
        Page 195
        Page 196
        Page 197
        Page 198
        Page 199
        Page 200
        Page 201
        Page 202
        Page 203
        Page 204
        Page 205
        Page 206
        Page 207
        Page 208
        Page 209
        Page 210
        Page 211
        Page 212
        Page 213
        Page 214
        Page 215
        Page 216
        Page 217
        Page 218
        Page 219
        Page 220
        Page 221
        Page 222
        Page 223
        Page 224
Full Text

Statistical Modeling and Simulation cf
Atmospheric Turbulence






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, co-chairman 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.



ACKNOWLEDGEMENTS................... ..... .............. ii

NOMENCLATURE.......................................... vi

ABSTRACT............... ............. ............... xi


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 Monin-Obukhov similarity
theory............................... 12

212 Atmospheric Diffusion...................... 14

2.2.1 Early Developments and the
K-theory............................. 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 Fokker-Planck equation........... 27

3.1.2 Stochastic differential equation..... 33

3.1.3 Comparison of the Fokker-Planck
equation with the diffusion equation. 37


3.2 The Growth Law ...........................

3,2.1 The Lagrangian similarity theory
for large times......................

3-2.2 The growth law for a Markov process...

THE STOCHASTIC MODEL............................

4,1 Description of the Model.....................


The mean time steps...........

The mean Lagrangian

. RESULTS..........................................

5.1 Computational Parameters....................

Estimation of Simulation Parameters.........

4,2.1 Neutral conditions.................... The mean velocity profile..... The variance of the velocities The variance of the
accelerations................. The mean time steps........... The mean Lagrangian

4,2,2 Non-neutral conditions................ The mean velocity profile..... The temperature profile....... The variance of the velocities The variance of the














5.2 Validation of the Model..................... 85

5.2.1 Neutral conditions.................... 86

5.2.2 Non-neutral 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 Non-Gaussian exponential function..... 153


6.1 Conclusions... ............................. 159

6.2 Recommendations... ......................... 161



B. COMPUTER PROGRAMS.................................... 170

LITERATURE CITED....................................... 206

BIOGRAPHICAL SKETCH...................................... 222


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

c concentration, ML3

c concentration distribution coefficient, ML-3
c ground level concentration at y=0, ML
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 y-direction

f deterministic function (Sec. 3.1)

f space correlation coefficient (Sec.

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 z-direction

g acceleration due to gravity, Lt-2

g random function (Sec. 3.1)

h deviation function from Gaussian distribution in the
z-direction without reflection

h' reflection term corresponding to h

H stack height, L









ml, mi2








Pl' P2


q1' q2


- deviation function from Gaussian distribution in the

- concentration integral, ML-2

- wave number, L1
2 -1
eddy diffusivity, L T-

- length scale of turbulence, L

- Monin-Obukhov 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
- 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, Mt-1.
- source strength, Mt

r generation of the species due to chemical reaction,

r exponentially distributed random variable with mean
of one

R gas constant (Sec., MI2t2T-1

R correlation coefficient

s normally distributed random variable with mean=0, and

t time

tE Eulerian microscale of turbulence, t

t sampling time, t

T temperature

TG gradient similarity temperature scale, T

T Monin-Obukhov temperature scale, T

T integral time scale of turbulence, t
u velocity in the x-direction, Lt1

S- friction velocity, Lt-1

Up convection velocity, Lt-

v velocity in the y-direction, Lt-

V- velocity vector, Lt1

V c volume of a cell L

w velocity in the z-direction, Lt-

x coordinate axis, downwind distance, L

X displacement of a particle in the x-direction, L

X displacement vector, L

y coordinate axis, cross-wind distance, L

Y displacement of a particle in the y-direction, L


z coordinate axis, vertical distance, L

z roughness length, L

Z displacement of a particle in the z-direction, L

Greek Letters:

a velocity profile exponent

a nt moment of particle displacement, Lnt-l

S- ratio of the Lagrangian and Eulerian integral time

S- constant in the log-linear velocity profile (Sec. 4.2.2)

8 exponent in Eq. (5.60)

Y exponent in Eq. (5.56)

r adiabatic lapse rate, TL-1

6 exponent in F equation (Sec. 5.4)
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
K von Karman constant

A mean value of the time steps
A integral length scale of turbulence, L

S- viscosity, ML- t-l
2 -1
v kinematic viscosity, L t-

dimensionless vertical distance (Sec. 5.4.2)

p density, ML-3

a standard deviation

T- lag time, t

T time scale of turbulence (Sec. 3.2), t

-1 -2
S- shear stress (Sec., ML t

-1 -2
T shear stress at the ground, ML t

dimensionless velocity scale

w velocity scale of particle motion, Lt-


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


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



Suresh B. Pahwa

March, 1975

Chairman: R.W. Fahien
Co-chairman: 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 Fokker-Planck equation. Comparison of the Fokker-

Planck equation with the diffusion equation showed that the deterministic

velocity used in the Fokker-Planck 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 cross-wind 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.


For non-neutral 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.



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 radio-active material

with sufficient accuracy and speed. Since the Clean Air Act of 1967,

some serious efforts have been made in the abatement of industrial


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 Fokker-Planck 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 Pasquill-Gifford charts for the neutral conditions so as to include


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.


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


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 three-dimensional continuum phenomenon, rotational, dissi-

pative, non-linear, 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

non-closed. 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


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 non-isotropic 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 k-5/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 non-linear 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


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 speed--where 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 one-dimensional

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 flux-gradient 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 velocity-velocity 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'.


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 Monin-Obukhov 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 (x-y) 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:

L = (2.8)


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)

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 turbulence--not 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 K-theory

The time averaged equation of continuity for a species in a


medium can be expressed in general as (Bird et al., 1960):

Dc 2
-+ V-V'c' = DV c + r (2.11)

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
velocity, and RL is the Lagrangian correlation coefficient. He assumed



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 K-theory. Velocity is assumed independent

of height and the eddy diffusivities in the y and the z-directions

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 z-directions 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 two-dimensional diffusion equation. Calder (1952) made a

contribution to the K-theory by expressing the eddy diffusivity as

K = K u, z (2.16)

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 K-theory 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
in the y-direction (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

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 Pasquill-Gifford

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

K = 1/3 4/3 (2.19)

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

non-neutral 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 well-known 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

three-dimensional grid, and the numerical solution may use a finite

difference scheme, a variational method or a particle-in-cell (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 Kirmse-Fahien 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 Kirmse-Fahien and the Bullin-Dukler

stochastic models is the representation of the eddy velocity. An

eddy velocity field continuous in shape is obtained from the

Kirmse-Fahien 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 non-zero 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 three-dimensional

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 Fokker-Planck



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 Fokker-Planck equation

Consider a continuous stochastic process X.. As we have seen

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

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


Since we are interested in p(X. ,t ), the probability density of the
1 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 n-1 11

P(Xil,t ; X. ,t ); .... X ,t ) (3.2)
1 2 n-1

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 n-1 n-2 1 n n-1


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 n-1
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 Chapman-Kolmogorov equation (Bhar-

ucha-Reid, 1960):

p(X,t Xi ,t ) = fP(Xi,tlX!,t')p(X!,t' IX ,to)dX' (3.4)
o o

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

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

to t in two stages.

Now let us consider an increment of time At. If the present

time in the Chapman-Kolmogorov 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


p(X ,t+AtIX. ,t ) =j P(X.,t+AtlXn-AX.,t)

p(Xi-AX.i,t X ,t )d(AX.)


Expanding the integrand on the right hand side of the equation in a

Taylor series in the x-direction about the point p(X.+AX.,t+At|X.,t)-

p(Xi,tJXi ,t ) we obtain:

p(X ,t+At iX.-AX.,t)p(X.-AX.,t|X. ,t ) =

= p(Xi+AXi,t+Atxi.,t)p(X.tIXi ,to) +

n! ) n


p(XitlXi ,to)]


Substituting Eq. (3.7) into Eq. (3.5), we obtain


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


[p(X +AX.,it+At Ixi.,t)p.

(Xi,tjXi ,to)]d(AXi)


[P(Xi L+!Xi't+AtlXi't)-


= p(X.,t Xi ,t ) / p(X.+AX.,t+AtlX,t)d(AX.) +

[P(Xi,t Xi ,to)
X [ o

1 1 1 1


From the definition of the probability density, we have


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

[p(Xi,tlX. ,to)
i0 o


Dividing Equation (3.11) by At and taking the limit as At 0:

At- 0

p(Xi,t+AtjXi ,t ) -p(X ,tXi ,t ) =
o 1 o

i (-1)n

8 n

[p(Xi,tli ,t )


p(Xi.+AX., t+At ., t)d(AX.)


At -+ 0

I n!


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

p(Xi,t+AtIXi ,to)-p(X.,tlXi ,t )
o o

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)


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

a. (Xi t) Lim -1 (AX) np(X +AXit+At X.,t)d(AX.)
n 1 At 1 1 1 1 1
1 At0


Substituting Equations (3.13) and (3.15) into equation (3.12), the

"forward Fokker-Planck equation" is obtained.

ap(xi,txi. ,to)


n Sai [ac (Xi,t)p(X.,ti x. ,t )]
n! a ni 1 n 0 0

S< X. < t

t > t



The boundary conditions are:

Lim p(X, tI X ,t ) = 6(X.-X )
t-0 o o

Lim P(X.,tlX. ,t ) = 0
Ixl 1 o

Equation (3.16) is also known as the forward Kolmogorov equation.

The three-dimensional Fokker-Planck equation is similar and the de-

rivation is identical to the one described above except for Equation

(3.6) which changes to a three-dimensional 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:

X. = X. + V At + a. m (3.18)
m m-1 m-1 m-1

where X. is the displacement of the particle at the end of the
(m-l)-- step. V. and a. are the velocity and the acceleration re-
1 1
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

entire past history of the particle. Therefore the displacement of

a particle in the Kirmse-Fahien 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-
lated as:

a. Li= A 1 f (AX.)p(X+AX.t+Atl.,t)d(AX.)
1. At-mo At f i

= E-- Xi.t = fi(Xi,t) (3.20)


a = (dX) Xt (3.21)
21 dt

The deterministic function f. is then the deterministic
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.


a (3.23)


S= a' (3.24)

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

= a = a + a' (3.25)

The Kirmse-Fahien 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
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 Fokker-Planck equation with the diffusion

The Fokker-Planck equation for a = 0 for n > 3 reduces to

S1 a2
S(Xti t) = (- p) + 2 (a p) (3.28)
at 0 3X. 1. 2 2 2.
o l l i

where X. is the particle displacement in the i-- direction. The
first moment in the i- direction a from Equations (3.20) and
(3.27) then is

= fi(Xit) = V. + f a. dt (3.29)

The first moment in the Kirmse-Fahien model a. is therefore
the deterministic velocity which is the sum of the average velocity

in the flow field and another term arising due to the mean Lagrangian


The physical meaning of a2. is easy to see from Equation
(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


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

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..
The Fokker-Planck 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)


al (Xit) = V. (xit) + f a (xi,t) dt (3.33)

Substituting Equations (3.32) and (3.33) into Equation (3.28) and

using 1/2 a2. = KL we have
1 1


ac (x ,t) -
('t) = [V.(xi.t) + a (x.,t) dt] c (xit)


+-2 [KL (x.,t) c(x.,t)] (3.34)
8x 1

The time-averaged diffusion equation, including the convec-

tive flow term, is

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-
sivity KL by equating Equation (3.34) and (3.35) and solving for
c a dt +a (L c)
1 1
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.
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
equation. If the diffusion process has Markovian properties, then

in general the moments of higher orders should be non-zero and should

be retained in the Fokker-Planck 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)

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.

dr a d (3.38)

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:

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-


2,(t) t [X(t)] -


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.

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 t-tv 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

d--t= bu (3.43)

dX t -) (3.
dt z

for t>> H

The concentration is assumed to be of the form

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




Based on the two above mentioned relationships, the exponent for

the plume width at the ground is found to be

m2 log -
z Z
0 0


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


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-


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


RL(T) = e (3.48)

where T is the integral Lagrangian time scale. If Z represents the

displacement of a particle in z-direction, 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

where a is the variance of the dispersion in the z-direction.
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


Substituting Equations (3.50) and (3.51) into Eq. (3.49), we have

K* z
( 0 -

dZ 2
2 = 2blu

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 2b K
z = [in --+i- 2
[lnZ +1)]2

a =bA 1+
z L AL l

Sx -
2 x AL
A 2 A- 1 + e L
L AL l

2 3
x x x
A1 +
L2 2!nZ 3!AL3

2 2b K
b 2
In( F + 1)1






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


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:

S2 3 -
b x x 2
S 1 ++(3.57)

S 12where


/ 2b K
b 1

ln(+ 1

z -
o H



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:



0 = x (3.61)
z 1
2 1 + x



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.
(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


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


The same form of the growth law is expected to be valid for

the dispersion in the y-direction as well. Thus

*a bx (3.65)
y H 1 + b x

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 3-1. Comparison of It curve calculated from the correlation coefficient data
with the Markov process curve fit for typical neutral conditions.


H 16




Markov process curve fit

0 II I I.
0 10 20 30 40 50 60 70
t, sec

Figure 3-2. Comparison of I curve calculated from the correlation coefficient data
with the Markov process curve fit for typical unstable conditions.






t, sec

Figure 3-3.

Comparison of It curve calculated from the correlation coefficient data
with the Markov process curve fit for typical stable conditions.


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 x-axis extends horizontally in the direction of the mean wind,

the y-axis is in the horizontal plane perpendicular to the x-axis

(cross wind), and the z-axis 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


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-
trained by integrating Eq. (3.25) twice in the i- direction over At .

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
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 m-1 m-1 m-l

The total velocity V. is the sum of the mean velocity V. and the
i i
m-1 m-1

fluctuating velocity. The model is capable of using any given form

of the mean velocity V. (x,y,z). The mean Lagrangian accelera-
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
in the computer to get the concentration values. The concentration

is calculated from the particle sampling according to the following


Q Ncell(ts)
c= Q cell (4.3)
p cell


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

It can be shown that

Q-- is the source strength represented by each particle.

-- is the mass of the pollutant represented by each sampling.

Qt Ncell
s cel is the mass of the pollutant given by all the samplings
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


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)


TE = RE (t) dt (4.7)

Expressing TL in terms of TE

TL =fR (st) d(St)

j RL (st) dt (4.8)

Using Eq. (4.5)

TL = RE (t) dt = E Tg (4.9)


EL \\



Using well known relations between the space and time macroscales

(Hinze, 1959)


AE = u TE

and therefore

A uT -
A T 6o
V V. v.
i L

where a is the intensity of turbulence in the i-- direction.




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)



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. 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


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

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-
tionality as -, we have

I u

du 1
dz K (z+z )18)

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



independent of meteorological conditions. The friction velocity

u* represents the magnitude of the wind velocity and is invariant

with height. The variance of the velocities

These are assumed to be equal to the Eulerian mean square


o2 = ( ) = (a= 2 (4.20)
v. v. v. E The variance of the accelerations

The Lagrangian acceleration is obtained from the Eulerian

acceleration by the Langrangian-Eulerian 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),

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-


In a later report, Cramer and Record (1969) were able to suc-

cessfully correlate all their energy spectra measurements by plotting
all of their data on a graph of dimensionless spectra 2 vs.
u* z
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:

e = 2v/ k2E(k) dk (4.23)

where k is the wave number.

2k n
k = (4.24)


E (k) = -- F(n)


The wave number k can be non-dimensionalized in a similar manner to

the frequency n if we let

* nz
n -


Multiplying Eq. (4.24) by z and substituting Eq. (4.26) into

it, we have

kz nz *
2 -- n
2f d


Similarly, dividing Eq. (4.25) by (u*z) and simplifying, gives

27TE u F *
=~ =~ F
2 2
u*z u*z


Substituting for E and k in terms of F and n in Eq. (4.23), we



0o + 2 2
E = 2v J 2n*z ) u z d ( 2



872r 2 2 ,2 ,
E 82 n F dn (4.30)

2 2
87r2 u* *
2 (4.31)


2 2
n F dn = (4.32)
0 87r v u

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


S= 2ru* B /2E* (4.33)
ai z vi 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

X = 0.5 (4.34)
t. a
1 a.


where the factor 0.5 is determined theoretically. The mean Lagrangian accelerations

The mean Lagrangian acceleration in the x-direction a is

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,

the mean Langrangian acceleration in the z-direction is also taken

as zero.

4.2.2 Non-neutral conditions The mean velocity profile

In the case of non-neutral conditions, one must use an addi-

tional length scale called the Monin-Obukhov 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-
editions, = 0 and therefore

p (0) = 1



The dimensionless shear i can be expanded in a Taylor ser-

ies about zero and for small values of /L (< 0.1), one can


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 Monin-Obukhov length is defined as follows in terms of

the turbulent heat flux:

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



where K and Kh are the momentum and heat turbulent diffusivities
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


Substituting Equations (4.41) and (4.42), Eq. (4.39) becomes

uL (Km dz


K u, (du dz)

T R idz



A gradient

similarity length LG is defined as

LKh dz
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

for zIL| < 0.1



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 log-linear 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

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


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)

where a is the slope ofa vs. ln( /LG) line and b is the intercept
at /LG = 1. For unstable conditions

a = a constant

= 0.124 (4.53) 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)

C dT dP = 0 (4.54)
p p






0.2 -

0 0


0 2
10-5 22
10 2





5 102

Figure 4-1. The velocity-profile exponent a vs. the stability parameter z** for stable

I r 1

Integrating Eq. (4.54) from the reference pressure PR to P, we have


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

Using the ideal gas law and neglecting the higher order terms, Eq.

(4.58) reduces to

0 = T + -- z (4.59)

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 temperature-mean 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

T dz u, (4.62)
G di K

The log-linear temperature profile equation then is

( z+z
e e= R6+T rin 0oo + z (4.63)


0 2 4 6 8

u, mps

Figure 4-2.

The potential temperature vs. the.mean velocity
for some Prairie Grass runs with stable conditions.







28 -

#152 0



22 I I I
0 2 4 6 8

u, mps
Figure 4-3. The potential temperature vs. the mean velocity for
some Prairie Grass runs with unstable conditions.

The gradient temperature scale TG is related to the Monin-Obukhov

temperature T* according to the following equation.

T = T
G Kh


Similarly the power law temperature profile is

T i
S= eR + I
R (0.1)t L G

(4.65) The variance of the velocities 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. 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


of the velocity derivative according to the following equations:


1 1
2 2

f(xi) = RE.


where f(x.) is the longitudinal space correlation coefficient. Also


2 i_2
vi x.


X. =0O

According to Taylor's hypothesis

t -u
3t ax



(u' x)2

= --2
--2 2
u tE.

Using Equations (4.14) and (4.21), we obtain

r2 a
ai BtE
1 E



 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). The mean .Lagrangian accelerations

The mean accelerations in the x and y directions are taken

as zero. The mean Langrangian acceleration in the z-direction is

obtained from the equation of motion in the z-direction.

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.


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


The deviation p can be calculated by expanding p into a two-dimen-

sional Taylor series about the point po as follows:

T+ P+ (4.76)
P Pp0 T Po


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

Substituting Eq. (4.81) into Eq. (4.72), we have

- 1 dP T
a = 1 + + gz
z p dz T P


Since the pressure and temperature profile are known one can there-

fore obtain an expression for a in terms of z. The temperature pro-
file for the neutral conditions is given by

T = 68 + r z
o R


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



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'^


S 0-- 1)


z> 10




Combining Equations (4.85) and (4.86) with Eq. (4.83), we have

z+z \
G n --G +
T o G




+ + **+
+R+ N z
R R 1

for z < L




( I + + 1 **)
RoD +~ R1

for z > LG


* z
z =-


** z
z L










N -

Similarly, the expression for the pressure deviation term becomes

P + dP
R R dz)

z P(- z
R dz


PR+ (dz)

J dP
P dz dz
S o

A linear approximation for the above equation is


z dP o 0]
PR dz dz




N N -
** 2 1 R
1+ R- z

d -


p op
0 0



Substituting Equations (4.57), (4.98) into Equation (4.82), we


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. (4-4). It can be seen that the variation in a
close to the ground is considerable.
close to the ground is considerable.



C 1 C **
N N- N z
R 2 2 iR

p S1+N **
0 1 R

R N2
1 -"- 2
Cp N


0 250 500 750 1000
z, m

Figure 4.4. The variation of different terms in the vertical mean Lagrangian acceleration
equation with height.


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, two-dimensional cell

structures were used at different downwind positions. Since perfect

symmetry is assumed in the y-direction, all the particle positions

are folded in the y-direction. 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 y-direction at the 50 m, 100 m and 200 m arcs, and in

the z-direction for two towers (y = 5 m and y = 18 m) are compared

in Figures (5-1) to (5-5), respectively, with the experimental re-

sults and with the simulation results of Bullin and Dukler (1974).

The concentration profiles, obtained from the Pasquill-Gifford equa-

tion (Turner, 1967) are also plotted on the same graph using the

standard deviations for the horizontal and vertical dispersion from

the Pasquill-Gifford 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

z-directions are plotted in Figs. (5-6) to (5-9). 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


O This work

-- Bullin & Dukler (1974) -
(class D)

0 0

- -I O


S .0



I 0@

Figure 5-1.

y, m
Concentration distribution in the y-direction for Project Prairie Grass
Run #45 at x = 50.



300 I-

200 r


1 i -.




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs