United States
Environmental Protection
Agency
Environmental Sciences Research
Laboratory
Research Triangle Park NC 27711
EPA-600 4-79-023
April 1979
Research and Development
A Lagrangian
Approach to
Modeling Air
Pollutant Dispersion
Development and
Testing in the
Vicinity of a Roadway
-------
RESEARCH REPORTING SERIES
Research reports of the Off ice of Research and Development, U.S. Environmental
Protection Agency, have been grouped into nine series. These nine broad cate-
gories were established to facilitate further development and application of en-
vironmental technology. Elimination of traditional grouping was consciously
planned to foster technology transfer and a maximum interface in related fields.
The nine series are:
1. Environmental Health Effects Research
2. Environmental Protection Technology
3. Ecological Research
4. Environmental Monitoring
5. Socioeconomic Environmental Studies
6. Scientific and Technical Assessment Reports (STAR)
7. Interagency Energy-Environment Research and Development
8. "Special" Reports
9. Miscellaneous Reports
This report has been assigned to the ENVIRONMENTAL MONITORING series.
This series describes research conducted to develop new or improved methods
and instrumentation for the identification and quantification of environmental
pollutants at the lowest conceivably significant concentrations, tt also includes
studies to determine the ambient concentrations of pollutants in the environment
and/or the variance of pollutants as a function of time or meteorological factors.
This document is available to the public through the National Technical Informa-
tion Service, Springfield, Virginia 22161.
-------
EPA-600/4-79-023
April 1979
A LAGRANGIAN APPROACH TO MODELING
AIR POLLUTANT DISPERSION .
Development and Testing
in the Vicinity of a Roadway
by
R. G. Lamb, H. Hogo, and L. E. Reid
Systems Applications, Inc.
950 Northgate Drive
San Rafael, California 94903
Contract No. 68-02-2733
Project Officer
R. E. Eskridge
Meteorology and Assessment Division
Environmental Sciences Research Laboratory
Research Triangle Park, North Carolina 27711
ENVIRONMENTAL SCIENCES RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
RESEARCH TRIANGLE PARK, NORTH CAROLINA 27711
-------
DISCLAIMER
This report has been reviewed by the Environmental Science Research
Laboratory, U.S. Environmental Protection Agency, and approved for publi-
cation. Approval does not signify that the contents necessarily reflect
the views and policies of the U.S. Environmental Protection Agency, nor
does mention of trade names or commercial products constitute endorsement
or recommendation for use.
-------
ABSTMCT
We have developed and demonstrated a microscale roadway dispersion
model based on Lagrangian diffusion theory. The model incorporates
similarity expressions for the mean wind and turbulence energy in the
atmospheric boundary layer, through which the effects of wind shear and
atmospheric stability are taken into account; a parameterization of vehicle
wake turbulence; and provisions for predicting the turbulence quantities
required to simulate second-order chemical reactions. Through simple
modifications, the model can be structured to treat particle settling,
deposition, and resuspension, as well as buoyancy of the effluent material.
Calm winds, winds parallel to the roadway, flows around depressed or elevated
roadways, shallow mixed layers, and transient or spatially variable meteor-
ological conditions can all be explicitly taken into account within the frame-
work of our modeling approach. The model requires minimum computer core
storage and fairly small execution times. It is not based on finite dif-
ference equations, and therefore, it is free of the truncation errors and
computational stability and convergence problems that models based on those
methods encounter.
We demonstrated the model by applying it to 52 of the 30-minute experi-
mental periods reported in the GM sulfate study. Of the 1040 predicted
values of the mean concentration of an inert material (SFg), we found half
of them to be within ±30 percent of the measured values. The overall
correlation coefficient was 0.91. The monitoring sites treated ranged in
distance from 2 to 100 meters from the edge of the roadway. The computer
time (but not core storage) required by our model is directly proportional
to the distance between the farthest receptor and the road. For the studies
reported here, the model required on the average 20 seconds of CPU time
on the CDC 7600 to simulate each of the 30-minute GM experiments.
m
-------
One finding from our tests of the model is that vehicle wake turbu-
lence is important in dispersing material near the roadway, especially
under stable atmospheric conditions. The great improvement in the accuracy
of the predicted concentrations achieved by adding just a simple para-
meterization of the wake indicates that concentration levels close to the
road are relatively insensitive to the intensity of ambient turbulence.
Further refinement and testing of this model is planned to improve
the wake turbulence parameterization, to evaluate in more detail the
accuracy of the method used to simulate ambient turbulence effects under
combined light wind and very unstable situations, to incorporate the non-
linear chemical reaction simulation capability, to make the simple modifica-
tions necessary to simulate fugitive dust dispersion, and to optimize the
computer code.
IV
-------
CONTENTS
Abstract i i i
Illustrations , vi
Tables viii
I. Introduction 1
II. Formulation of the Model 7
A. General Formulation 7
B. Determining the Function J and Its Time Integral 13
C. Application to Photochemical Pollutants 18
III. Model Testing 21
IV. Meteorological Inputs to the Model 27
V. Results of Preliminary Model Validation Experiments 52
A. Simulations That Neglect Wake Effects 52
B. Simulations That Include a Simple Parameterization
of Wake Effects 55
VI. Summary 68
Appendices
A. A Scheme for Simulating Particle Pair Motions
in Turbulent Fluid 69
B. User's Guide to the Lagrangian Highway Model 91
References 104
-------
ILLUSTRATIONS
Number page
1 Schematic Illustration of the Relationship Between the Lagrangian
and Eulerian Approaches to Turbulent Diffusion 4
2 Coordinate System Used in the Roadway Model To Compute Concentra-
tions at an Arbitrary Receptor Site 8
3 Schematic Illustration of Some of the Parameters Required and
the Method Used TO Calculate p 15
V
4 Example of the Dependence of J^ on the Number of Realizations K. . . 17
_ p
5 Simple Particle Trajectories with u = 3.0, u1 =1, and
w7? = 0.6 (MKS Units) 22
6 The Root-Mean-Square Particle Dispersion Component o as a
Function of Travel Time for Two Values of the Integral
Time Scale T 23
7 Comparison of Roadway Model Predicted Cross-Wind Integrated
Concentration with Gaussian Plume Formula Predictions at a
Distance of 100 m Downwind of a Point Source at 50 m Elevation ... 24
8 Comparison of Roadway Model Predicted Cross-Uind Integrated
Concentrations with Gaussian Plume Formula Predictions at a
Distance of 100 m Downwind of a Point Source at 5 m Elevation. ... 25
9 Coordinate System Used To Express Mean Wind and Turbulence
Velocities 29
10 Observed Versus Predicted Mean Wind Speed 36
—x-1/2 -VI / 2
11 Observed Versus Predicted Values of u'^ and wl<: 39
-^1/2 —?l/2
12 Observed Versus Predicted Values of u' = u1^ and w1 = w1^ ... 42
13 Observed Versus Predicted Values of u' and w' 44
14 Plots of the Ratio of the Wind Speed at 10.5 and 1.5
Meters and the Drag Coefficient u*/u(10.5) Given by
Similarity Theory as Functions of 1/L 47
-------
Number Page
15 Comparison of the u* and 1/L Values of Table 1 with Those
Derived Using the WSR Method 49
16 Comparison of Observed SF5 Concentrations with Those
Simulated by the Roadway Model with Vehicle Wake Turbulence
Neglected: All Experiments Included 54
17 Comparison of Observed SF5 Concentrations with Those
Simulated by the Roadway Model with Vehicle Wake Turbulence
Neglected: Data Points for Day 295 (H < 1 m/sec, Ri > 0.25)
Omitted 56
18 Comparison of Observed SFs Concentrations with Those
Predicted by the Roadway Model with Wake Turbulence Effects
Parameterized: All Experiments (1100 Data Points) Included . . 59
19 Comparison of Observed SF5 Concentrations with Those
Predicted by the Roadway Model with Wake Turbulence Effects
Parameterized: Results Displayed for Stably Stratified
Flows Only (580 Data Points) 61
20 Comparison of Observed SFs Concentrations with Those
Predicted by the Roadway Model with Wake Turbulence Effects
Parameterized: Results Displayed for Unstably Stratified
Flows Only (520 Data Points) 62
21 Comparison of Observed Peak SFg Concentrations for Each of
the 55 GM Experiments with the Values Predicted by the
Model (Wake Effects Included) 64
22 Distribution of Fractional Error 67
2 1/2
A-l Stceamwise Component U2 - £0)x of the rms Separation of
Simulated Particle Motions Under the Conditions u - 3.0,
u7"2 = 1, and w7? = 0.6 for Two Values of the Parameter aQ ... 89
VII
-------
TABLES
Number Page
1 Estimated Values of u* and_L for the GM Experiment Periods
(zn - 5 cm) and the Angle 6 of the Mean Wind Relative
to the Road 33
2 Model Performance Statistics 66
B-l Format for Card 1 93
B-2 Format for Cards 3 and 4 94
B-3 Format for the Monitor Site Cards 96
B-4 Format for the First Meteorological Data Card 96
B-5 Format for the Second Meteorological Data Card 98
B-6 Format for the Source Location Card 98
B-7 Format for the Miscellaneous Data Cards 99
B-8 Sample Input Deck Used in the GM Simulations 100
vm
-------
I INTRODUCTION
The General Motors (GM) Sulfate Dispersion Experiment conducted in
late 1975 and reported by Cadle et al. in 1976 revived interest in roadway
diffusion modeling. The study not only raised a controversy about the
adequacy of existing roadway models, but also collected a unique set of
turbulence and concentration data in a roadway environment under well-
monitored conditions. The work described in this report is in part an
indirect outgrowth of the GM study: Although the basic theory of the
model described here was developed several years ago, it was implemented
and tested under contract to EPA as part of the roadway model development
activity at EPA that the GM study encouraged.
Prior to 1976, the most widely used model for performing environmental
impact assessments of roadways was HIWAY, a model based on the Gaussian
plume formula (Zimmerman and Thompson, 1975). Comparisons made by Chock
(1977) of the predictions of HIWAY and concentration values measured during
the GM study revealed that HIWAY seriously overpredicted concentrations,
especially during stable, light wind conditions. The generally poor per-
formance of HIWAY was not surprising, considering the host of simplifying
assumptions inherent in the Gaussian plume formula on which it is based.
In addition to neglect of streamwise diffusion, an assumption that is
invalid during light winds, other simplifying assumptions include neglect
of vehicle wake turbulence, dismissal of wind shear effects, and use of
interpolation to obtain values for the dispersion parameters a and GZ
from some arbitrarily chosen initial value at the source and the values
measured approximately 100 meters downwind. Furthermore, the plume
formula is not applicable under nonsteady or spatially inhomogeneous
situations, it cannot describe dispersion that occurs around depressed
or elevated roads, and it is limited to the treatment of inert pollutants.
-------
Seeking to develop roadway models free of the restrictions inherent
in the Gaussian plume formula, many investigators have turned to the dif-
fusion equation. Among these are Danard (1972), Ragland and Pierce (1973),
Egan, Epstein, and Keefe (1973), and Eskridge and Demerjian (1977).
However, the diffusion equation is not a rigorously sound basis for multi-
point or line source models. In deriving the K-theory diffusion equation
from basic principles, one must assume that the smallest length scale of
the spatial variation in the mean concentration distribution is much
larger than the Lagrangian integral length scale of the turbulence (see,
for example, Lamb, 1973).
The effects of this assumption and the limitation it imposes can be
seen clearly in the special case of stationary, homogeneous turbulence.
In this case, the eddy diffusivity K is a constant and the diffusion
equation can be solved analytically. Its solution for a point source has
2
the standard Gaussian form with a = 2Kt, where t = x/u, x is distance
downwind from the source, and u is the mean (uniform) wind speed. This
expression is consistent with Taylor's (1921) classical result for diffu-
sion at large distances from a source, but it is at odds with the observed,
theoretically derivable relationship o = u^ t2 that applies when t is
small. Here, u"^ is the variance of the turbulent velocity fluctuations.
Thus, the assumption of smooth concentration fields invoked in the deriva-
tion of the diffusion equation implicitly filters out the turbulence
processes that are important during the initial period of diffusion from
2 2
a point source and that are responsible for the o ^ t relationship cited
above. By assuming that the diffusivity K is a function of travel time t,
viz., K = 1/2 u1^ t, one can make the diffusion equation compatible with
the observed behavior of point source dispersion. But this modification
can be used only if a single point or line source is present. Otherwise,
no unique value of K(t) or K(x) exists.
Teske and Lewellen (1978) have developed a roadway model using a
second-order turbulence closure technique. It is not clear to what extent
this approach circumvents the problem just cited. Two disadvantages of
-------
diffusion models based on second-order turbulence closure is that they are
complicated and costly to operate. Also, they suffer from several problems
that any model based on finite difference methods incurs. One of these is
that when the winds are very light or are blowing parallel to the road, the
pseudo-diffusion associated with many finite difference forms of the dif-
fusion equation can spread material over much larger areas than ambient
and vehicle wake turbulence together are capable of doing. Finite dif-
ference schemes that minimize pseudo-diffusion do so at increased costs
in computer time and storage requirements. A second problem with finite
difference models is that the number of grid points required, and hence
the computer time needed, increases with the square of the width of the
region treated. This occurs because the top of the model tr.ust be at
least several times the maximum value of o for the area of interest,
and as we pointed out above, a ^ x initially. As a result, it is
difficult to achieve the fine spatial resolution needed around the
roadway without making the core and machine time requirement of the model
excessive. Also, it is awkward to apply finite difference methods to
areas having complex shapes, such as those of elevated and depressed
roadways. Finally, one must always bear in mind that the finite dif-
ference equation is a model of the differential equation it represents
and, unless care is taken, under some conditions the finite difference
equation might possess solutions that are radically different from any of
those of the diffusion equation.
All of the problems discussed above can be circumvented by casting
the roadway model in Lagrangian form. Figure 1, which is taken from Lamb
(1978a), illustrates the way in which this approach is related to the
Gaussian plume and K-theory models discussed above. The theoretical
aspects of the model are developed in Chapter II. We discuss below sev-
eral important features of the Lagrangian approach that make it ideally
suited to simulating diffusion and even some chemical reactions over
short distances (up to several hundred meters) from point or line sources.
-------
FIGURE 1. SCHEMATIC ILLUSTRATION OF THE RELATIONSHIP BETWEEN THE
LAGRANGIAN AND EULERIAN APPROACHES TO TURBULENT DIFFUSION
-------
First, Lagrangian models require no grid network; as few or as many
receptor points as desired can be used, and these can be arbitrarily dis-
tributed in any configuration over the area of interest. The computer
core needs are therefore minimal. Also, the absence of a grid network
and of finite differencing schemes makes the computational task of modeling
dispersion around elevated or depressed roadways (or even that of modeling
transport and diffusion in complex terrain) relatively simple.
Secondly, the Lagrangian approach is essentially free of the restrict-
ing assumptions (mentioned earlier) that hinder the plume and K-theory
models in certain instances. The roadway model we develop here can
explicitly account for the following effects: wind shear; particle settling,
deposition, and resuspension; buoyancy of the effluent; calm winds or
flows parallel to the roadway; vehicle wake turbulence; atmospheric
stability; and time or space variability in the meteorological or source
emissions conditions. The temporal evolution of the dispersion from the
2 2
3 ^ t1" regime to that of o ^ t described above is also properly duplicated.
Thirdly, when used in conjunction with the closure approximation devel-
oped by Lamb and Shu (1978), the Lagrangian roadway model can treat pollu-
tant species such as NO, NOp, 0~, and others, taking into account the
effects on chemical reaction rates of turbulence-induced concentration
fluctuations.
Finally, this particular modeling approach utilizes readily measurable
Eulerian turbulence statistics, primarily the variances of the velocity
fluctuations, rather than Lagrangian statistics like the a parameter used
in puff and plume models. In this report, we show that available empirical-
theoretical expressions for the vertical profiles of mean wind and the
required velocity variances are accurate enough for obtaining good concen-
tration predictions through this modeling technique using only the fric-
tion velocity u*, the Monin-Obukov length L, and an estimate of the
Lagrangian integral time scale of the turbulence as input parameters.
-------
The remainder of this report discusses the formulation of the model
(Chapter II), testing of the model through comparisons of its prediction
with those of another model (Chapter II), the meteorological inputs to
the model, and the results of our preliminary model validation experiments.
Appendix A describes a scheme for simulating particle pair motions in
turbulent fluid that is the heart of the roadway model. Appendix B
describes the input variables required by the model.
-------
II FORMULATION OF THE MODEL
This chapter discusses the theoretical formulation of the roadway
model and its application to photochemical pollutants.
A. GENERAL FORMULATION
For simplicity, the roadway model is cast in a two-dimensional form.
The process by which this is achieved is illustrated in Figure 2. For a
roadway link of sufficient length (to be established later), oriented as
shown in Figure 2, the mean concentration at the receptor shown in the
figure is given by:
t »
= ( f p(x,0,z,t|xs,y',zs,tl)S(y',tl) dy' dt'
(1)
where z is the effective height of the vehicle emissions; x denotes the
location of the road in the chosen frame (only a single lane of traffic
is considered here); S is the function describing the distribution of
vehicle emissions; and p(r,t|rs,t') is the probability density that a
particle released at r1 at time t1 will be found at the point r at time t.
In Lamb (1978a), Eq. (1) is derived from basic Lagrangian principles. Note
that Eq. (1) is the generalized form of the solution of the inhomogeneous,
ensemble-averaged mass continuity equation of an inert substance, and p is
the Green's function of that equation. For traffic of uniform speed v$,
we have:
N
s(y'.t') = £ Qn6(y' -y0n - vnf) , (2)
where Qn is the emissions rate (mass/time) of the n-th vehiclo; N is the
total number of vehicles operating on the link; yQ is a random variable
representing the position at t = 0 of the n-th vehicle; vp = vg if the
-------
n-th vehicle is moving toward +y, vn = -v otherwise; and 6 is the
delta function.
RECEPTOR
COORDINATES = (x,0,z)
H = mean wind vector
u = component of V,, nonnal to roadway
v = component of VH parallel to roadway
0 = coordinate origin; y axis parallel to road
with receptor lying in x-z plane
FIGURE 2. COORDINATE SYSTEM USED IN THE ROADWAY MODEL TO COMPUTE
CONCENTRATIONS AT AN ARBITRARY RECEPTOR SITE
Owing to the symmetry of the roadway about the x-z plane, we can
assume that the component of diffusion parallel to the roadway is negli-
gible, in which case,
p(x,0,z,t|x ,y',z ,t') = p(x,z,t|xs,zs,t')6[y' + (t - t')v] , (3)
-------
where v Is the component of the mean wind parallel to the road (see Figure
2). Substituting Eqs. (2) and (3) into Eq. (1) and evaluating the y' inte-
gral , we obtain
- Z Qny P(x,z,t|xs,zs,t'h[t'(v - vn) - tv - y0n] df
(4)
Let a = (v - vn)t' = Vnt'. Then Eq. (4) can be written in the form:
N _V_t
- vt -
,o 'I - •**»„/"
r n
= S Qn / p(x,z,t|xs,zs,^-
n=l * ^
where
-^n) , 1f 0 0, or if
0 > vt + ynn > V t
^un n
and Vn < 0;
0 , otherwise. (5b)
To arrive at Eq. (5a), we made use of the following well-known property of
the delta function:
/*b
/ f(x)6(x - xQ)dx = f(x0) , a < XQ < b
*a
It is perhaps instructive to comment briefly on the physical meaning
of Eq. (5). The delta function in this equation is the effective source
strength function. Keeping in mind that 6(5) is nonzero only when ? = 0
and that yn is a random sequence describing vehicle locations at t = 0,
-------
we see that the source function S in our roadway model is a random sequence
of instantaneous point sources distributed in time in a manner determined by
the combination of vehicle and wind speed components v and v. For example,
if v = 0, i.e., if the winds are either calm or blowing perpendicular to
the road, then the delta function in Eq. (5) is nonzero only when
t'vn + yQn - 0. By definition of yQ , this occurs only at the instant
t1 when the n-th vehicle passes through the x-z plane of the receptor.
This obvious result occurs because, when v = 0 and the roadway is symmetric
about the receptor plane, only those pollutants released in that plane can,
in effect, influence that receptor. As a second example, suppose that all
vehicles are stopped, v = 0, but v i- 0. In this case, the delta function
in Eq. (5) is nonzero, for a given t, when (t - tl)v = -yfi . In other
words, emissions released at time t' from the n-th vehicle affect the
receptor at the tiir
the receptor plane.
receptor at the time t = t1 - (yn /v) when those emissions pass through
Consider the special case where V = v - v =0, which represents the
situation where the speed and direction of the n-th vehicle are identical
to those of the wind component parallel to the roadway. In this case, the
equation is best obtained from Eq. (4) rather than from the transformed
Eq. (5). We then have:
N J:
<(c(x,0,z,t))= £ Qn6(v t + y ) / p(x,z,t|x ,zs,t') dtl . (6)
n=l JQ
This is simply the superposition of N two-dimensional plumes (recall that we
are neglecting diffusion in the y-direction), each located at the y-position
of its source at time t, i.e., at y = y~ + vnt. The n-th plume passes
through the receptor plane y = 0 at time t = ~yQn/vn- Tnis relationship
reveals the obvious result that, if the vehicles are northbound, i.e.,
v (= v) > 0, then only those whose initial locations are south of the
receptor plane, i.e., yn < 0, will affect concentrations at the monitoring
site.
In most problems of interest, it is the time-averaged concentration
near the roadway that is of concern. The averaging time T should be long
10
-------
compared with the time scales of mixing and transport from the road to the
farthest monitor, but short, compared with the smallest time scale of tem-
poral changes in the vehicle flow rate and meteorology; otherwise the time-
averaged concentration would depend on the length of the averaging interval
A time average is denoted by an overbar ( ) and defined as:
t+(T/2)
(7)
on t is a consequence of the dependence of c on vehicle
flow rate and meteorology, which are functions of the time of day, t. From
here on, we drop the reference to t in
-------
where n is the flow rate (vehicles/time) of vehicles of type m, all of
which have the same emission rate Q , speed v . and location (x ,z ) in
m rn sm sm
the receptor plane. For example, automobiles in one lane might represent
the group designated by m = 1, trucks in that lane might constitute m = 2,
autos in another lane might be represented by m = 3, and so on. [To obtain
Eq. (9), note that the number N of vehicles that pass through the receptor
plane in a period T is N =X)fLi ii T.]
Equation (9) is identical to that of the mean concentration produced by M
steady, continuous line sources of strengths n Q (|v - v |)~ . The depen-
dence of the effective line source strength on the wind speed component v
is noteworthy because it has not been taken into account in previous line
source roadway models. The error produced by this omission is small in
cases of light winds where |v | » |v| ; but on congested roadways or in
situations where there is a strong component of the wind parallel to the
road, errors of 20 percent or larger can result in the estimates if v
is neglected in the source strength calculation.
Looking back over the steps that led to Eq. (5), we find that this
equation becomes invalid if a significant fraction of the emissions that
affect the receptor arise from portions of the roadway that are not parallel
to the y-axis. A condition on the length of roadway that must be straight
in the vicinity of the receptor to render Eq. (5) accurate can be derived as
follows. Let
Kc
(.0,
where (XD,ZD) denotes the receptor coordinates,
K K
uc = ^(7+ax) , Ola)
wc=^(z+oz) , (lib)
and (x,z) is the centroid of an instantaneous cloud of emissions and o and
/\
o are (roughly) the cloud dimensions. [See, for example, Monin and Yaglom
12
-------
(1971) for similarity theory expressions for x, z, o , and o .] Under these
A £-
definitions, T. is a measure of the time required for material released at
(x . z) to reach the receptor (XD, ZD) through the combined actions of
S S K K
transport and diffusion. Thus, the roadway in the vicinity of the receptor
must be straight for a distance L that satisfies the condition
Ls » vTL , (12)
because otherwise significant quantities of material could reach the receptor
from points whose projection on the receptor plane are different from the
assumed point (x , z ). Subject to this constraint, Eq. (8) can be used as
the basis of a roadway model in general situations, and Eq. (9) can be used
in cases in which the vehicle and wind speed components are such that the
line source approximation is valid [see the discussion preceding Eq. (9)].
B. DETERMINING THE FUNCTION J AND ITS TIME INTEGRAL
The kernel p that enters in the definition of J [see Eq. (5b)] is the
probability density that a material particle released at the point Ug.O
at a specified time, say t1, will be at the point (x,z) at the later time t.
In situations in which the meteorology is quasi-steady, as we have already
assumed, the turbulence can be considered to be stationary, in which case
the kernel p is a function only of the travel time T = t - t1 and not of t
and t1 jointly. In this case, by virtue of the definition of p, the quantity
P(X,Z,T|X ,z ,0) dx dz is equal to the fraction of all realizations of the
turbulence in which a particle released at (x ,z ) is located after a travel
time T in the box of dimensions (dx.dz) centered at (x,z). That is,
ry 7 iy 7 n\ - lim 1 V V(x,z.T|x ,z ) , (13a)
,X,Z,T x ,z ,u; - . ^ K s s
S S N->-«> N ax aZ i._i
where
, if a particle released at (xs,zs)
is in the box (dx,dz) centered at
I , t (x,z) after a travel time T on
*kvx,z,T|xs,zsJ = <^ realization k;
, otherwise.
13
-------
Combining Eqs. (13a) and (5b) and integrating with respect to T, we
obtain:
T ]/
/ lim 1 A
J(x,z,T|x ,z ) dt = K™ i/ Hy H7 2^ At,(x,z|x ,z ;T) , (14)
S S IV'33 i\ UA Ui. . _-, K. b b
u
where At,( x.z! x ,z ;T) is the total time spent during the period T by a
K
particle i
at (x,z).
particle released at (x ,z ) on realization k in the box (dx,dz) centered
Equation (14) provides an operational procedure for implementing the
roadway model numerically: If the random turbulence velocity fields could
be simulated by a Monte Carlo type of process, then particles could be
released at any desired locations (x ,z ) in the simulated flows, their
motions could be tracked, and Eq. (14) could be evaluated computationally.
The resulting values could then be used in either Eq. (8) or Eq. (9) to
compute mean concentrations at any desired locations.
We can illustrate the method in more detail with the aid of Figure 3.
For the given roadway site and meteorological conditions, we specify the
profiles of the mean wind u"(z)--note that only the component u normal to
the roadway affects the kernel p, the profiles of the variances u'2(z) and
—P
w (z) of the ambient turbulent velocity components in the receptor plane,
and the profile of the Lagrangian integral time scale r(z). The functional
forms used for each of these profiles under various meteorological conditions
are listed later in Chapter IV. To account for the effects of vehicle wake
turbulence, it is also necessary to specify the speed and drag coefficient
of the vehicles and certain other parameters descriptive of their wakes. At
present, there are no accurate theoretical or empirical descriptions of the
wake turbulence generated by automobiles moving along a roadway. Consequently,
only approximations of the wake effects can be incorporated into the disper-
sion model. The method we have chosen is to express the variances of the
lateral and vertical components of the wake turbulence velocities as functions
14
-------
of the energy dissipated by drag forces on the vehicle and distance from
the vehicle and then to estimate the integral time scale T .of the wake
turbulence. The last quantity and the prescribed wake velocity variances
are then used to simulate particle dispersion by the wake in the same
manner (described later) in which the velocity and time scale properties
of the ambient turbulence depicted in Figure 3 are used to describe dis-
persion by boundary layer turbulence. The technique is described below.
RECEPTOR
SAMPLE VOLUMES
SIMULATED PARTICLE
TRAJECTORY
FIGURE 3. SCHEMATIC ILLUSTRATION OF SOME OF THE PARAMETERS REQUIRED
AND THE METHOD USED TO CALCULATE p
Let (x ,z ) be the release point at which we wish to evaluate the
j o
density function p. A particle is released from this point at time T = 0.
In the atmospl
are given by:
In the atmosphere, it acquires velocity components u (T) and w (T) that
(15a)
W (T) =
= W[ZP(T).T] + W"(T)
(15b)
where (x , z ) denotes the projection of the particle's coordinates on the
receptor plane; tl is the mean wind speed component normal to the roadway
15
-------
(depicted in Figure 3); u1 and w1 are random velocities representing the
ambient turbulence; and u" and w" are random velocities representing the
vehicle wake turbulence. If there is a nonzero mean vertical wind component
w or if the particle is positively or negatively buoyant, Eq. (15b) must
be modified accordingly.
In the numerical simulation, the random velocity components u1 and w1
are generated at each time step by an algorithm, described in detail in
? ?
Appendix A, that utilizes the variance profiles u' and w' and the
Lagrangian time scale T of the ambient turbulence. The random wake velo-
cities u" and w" are generated by the same process, except in this case the
energy and Lagrangian time scales of the wake turbulence are used. (Further
details are given in Chapter V.) Once values for all the random velocity
components have been specified, these are used together with the given mean
wind profile u in Eq. (15) to calculate the particle trajectory
[x (T),Z (T)] for 0 <. T <. T. As mentioned earlier, the integrating time
period T must be large enough for any particle to pass just beyond the
location of the most remote monitor.
At those sites at which the predicted particle trajectory passes
through the sample volume (dx,dz) of a monitor, the time At, that the
particle spends in that volume is determined. Here the subscript 1 signifies
that this is the first realization of the turbulence processes.
All of the steps outlined above are now repeated, starting with the
release of a new particle at (x ,z ) and using a new set of random velocity
j o
components u1, w1, u", and w" that are uncorrelated with the first set. From
the new realization of the particle trajectory, we determine the residence
time At2 in each receptor's sample volume. By continuing this process K
times, we obtain the sequence of At. values required in Eq. (14) to evaluate
the time integral of 0 that enters into the equation for the mean concentra-
tion. To determine how many realizations K are needed to render a finite
sum (over K realizations) an acceptable approximation of the infinite sum
of Eq. (14), we use a method in which intermediate approximations of the J
integral are stored and analyzed for trends that reveal convergence.
16
-------
Specifically, let
K
Atk
represent the estimate of the J integral given by K realizations of the
_K
particle trajectories. Figure 4 shows how 0 might vary with K. The
estimate usually varies rapidly as K increases from zero, but it becomes
_i/
increasingly steady as K becomes large and J converges to the true value
K
of the J integral. As a measure of the variability of J at realization
K, we use the standard deviation OJ(K) of the n + 1 values J*, 0K"AK, ...,
LTn , where n is any integer greater than 2 (we use 5) and AK is some
sampling interval in K (we use AK = 100). We consider that CT is an accept-
i/
able approximation of the true value of the J integral when oj(K) <. AiT,
where X is some fraction, e.g., A = 0.01. In our studies with n = 5 and
I/
AK = 100, we have found that the convergence condition a,(K) <. \j is
j
satisfied within K = 500 realizations at points near the core of a point
source plume, while K = 2500 may be necessary to satisfy the criteria at
I/
the edge of the plume where J is only one-tenth the centerline value or
less. We have also found that a great reduction in the computer time
required by the model can be achieved with virtually no loss of accuracy
by setting A = 0.5 rather than 0.01. (We were surprised to find that by
using A = 1.0, the overall performance of the model simulations of the GM
experimental data deteriorated by only 15 percent, as measured by the
change in the rms error, while the computer time required dropped to 1/30
of that necessary to run the model using \ = 0.05.)
K
The value of J derived from carrying out this procedure for each
monitoring site is used in either Eq. (8) or Eq. (9) for the J integral
to obtain the mean concentration. Chapter III discusses the agreement
between the prediction of this model and those of the Gaussian plume
formula for the meteorological conditions under which the latter is
known to apply.
17
-------
K'-n-K
FIGURE 4. EXAMPLE OF THE DEPENDENCE OF T ON THE NUMBER OF
REALIZATIONS K. The difference between 3* and the
time value of the time integral of J is estimated
from (n + 1) moving of U^ as described in the text.
APPLICATION TO PHOTOCHEMICAL POLLUTANTS
The roadway model can be applied to photochemical pollutants by
using the closure technique developed in Lamb and Shu (1978). To facili-
tate this operation, we assume that pollutants are either pseudo-inert
or infinitely fast reacting. This division is prompted by the fact that,
in most cases of interest, T is of the order of several minutes, which is
a long time compared with the reaction time of species such as NO, 03,
and sulfates around a roadway, but a short time compared with other
species such as hydrocarbons. The category a given species falls into
must be decided case by case, but we anticipate that around roadways, NO,
03, and sulfates can generally be treated as infinitely fast reacting.
In this case, the equations governing the mean concentrations of these
species are (from Lamb and Shu, 1978):
cNO(r,t) = max JO , [cINO(r,t) -
= X0 (r,t) - cINQ(r,t) + cMn(r,t)
m^'~
(17)
(18)
(19)
18
-------
Equation (19) assumes that sulfate is emitted by motor vehicles rather than
formed in the atmosphere through the oxidation of SC^. Here XQ is the
ambient ozone concentration, assumed to be given; cfjNO is the mean
concentration of NO in the absence of any chemical reactions [C\NQ can be
calculated direcely from Eq. (8) or Eq. (9)]; and £NQ is a mixing parameter,
which we approximate by (see Table 1 of Lamb and Shu, 1978):
_2
r (r t) - CINO(!>t) (20)
C^'I; -— ' { '
In Eqs. (17) through (20), r = (x,z) denotes the receptor location.
In summary, Eq. (8) is the roadway model equation for the mean
concentrations of all pseudo-inert pollutants. This same equation is
also used to compute quantities like CJNO that enter in Eqs. (17)
through (20). The model equations for NO, 0.,, and sulfate, which we
consider to be rapid reactants, are Eqs. (17) and (19). Similar equations
can be written for other rapid reactants. The evaluation of Eqs. (17)
and (19) requires that the mixing parameter c, be determined.
The Lagrangian theory from which Eq. (8) was derived (Lamb, 1978a) also
provides an expression for c2 that can be used to evaluate ;. We have:
_ N N t t
c (r,t) = Q /]/J I I p«(x,z,t;x,z,t x ,z ,t';x ,z ,t")
~" M^^MP ^^^^^> V • £
n=l m=l 0 0
*[<$ t'(v - vsn) - tv - y0nj-6[t"(v - vsm) - tv - y0 1 dt1 dt"
(21)
where p« is the joint probability density for the displacements of a pair
of particles. Note that p is a marginal density of P2 in that the former
can be obtained by integrating p2 with respect to the position coordinates
of one particle over all space. The method described earlier for
19
-------
determining the J integral can be used to obtain the integrals over p2 in
Eq. (21). However, in this case we must release pairs of particles from
the source point (x ,z ) and determine the total times that both particles
-> 5
are simultaneously within the sample volumes of each receptor. The algo-
rithm described in Appendix A that we use to generate the random velocity
components u1, w', u", and w" of single particles can also be used to
simulate the velocities of pairs of particles where motions are correlated.
In this role, the algorithm requires (in addition to the turbulence velo-
city variance profiles and Lagrangian integral time scale) the energy
dissipation rate profile e(z) of the ambient turbulence and comparable
information for the vehicle wakes. Some results of particle pair simula-
tions and mean square concentration prediction are presented in
Chapter III, but applications of the roadway model to photochemical pollu-
tants were not treated further in this study.
20
-------
Ill MODEL TESTING
Since the model described in the previous chapter is somewhat new,
it is appropriate to compare its predictions with those of other diffusion
models in current use. In particular, by comparing the predictions of our
model with concentration estimates obtained from the Gaussian plume formula
for situations in which this formula is known to apply (namely, under
conditions with no wind shear, steady meteorology and source emission
rate, no chemical conversion or deposition, and so on), we can verify our
approach as being basically sound. This chapter discusses these verifica-
tion efforts.
We set the model up to simulate dispersion from a point source in a
uniform flow (no shear) of speed u = 3.0 m/sec over a flat surface. The
—5" 2 2
variances of the ambient turbulence velocities were set to u = 1 m /sec
—5" 2 2
and w1' = 0.61 m /sec . Figure 5 shows 20 simulated particle trajectories
from a source 50 m high for each of two values of the time scale T.
Note that as T •*• 0, the value representative of Brownian motion, the
trajectories acquire the erratic character associated with molecular scale
movement. The root-mean-square (rms) particle dispersion component o
X
corresponding to the cases presented in Figures 5(a) and 5(b) are plotted
112
in Figure 6. It is evident that the a ^ t and a ^ t ' regimes predicted for
short and long travel times t, respectively, by Taylor (1921) are reproduced
by the scheme we use to generate the turbulence velocities. The pre-
dicted profiles of mean, cross-wind integrated concentration under the
same conditions such as those shown in Figures 5(a) and 5(b)
are shown in Figure 7 along with the corresponding values given by the
Gaussian plume formula. A similar comparison is given in Figure 8 for
a source height of 5 m. The excellent agreement between the model and
21
-------
cO
o.
CSI
100
2DO
""t" 1 ^T I I I I"
3DO
I 1 I I I I I
4oa>
o
rvi
I §4-
o
4-o
200
X (meters)
(a) T = 3 Seconds
zoo
X (meters)
(b) 31 = 30 Seconds
FIGURE 5. SAMPLE PARTICLE TRAJECTORIES WITH u" = 3.0,
u7? = 1, AND w7? = 0.6 (KKS UNITS)
22
-------
100
9
8
7
6
C
o
o
0)
00
X W
3 »
8
7
6
21 = 3 sec
SLOPE = 1/2
1 I I I I I I
7 3 4 i 6 7 1
7891 2 3 4S67«910
t (seconds)
9 WO
FIGURE 6. THE ROOT-MEAN-SQUARE PARTICLE DISPERSION COMPONENT a
AS A FUNCTION OF TRAVEL TIME FOR TWO VALUES OF THE
INTEGRAL TIME SCALE T
23
-------
90
80
70
60
50
'40
30
20
MODEL PREDICTION
CALCULATED VALUES FROM
GAUSSIAN PLUME FORMULA:
h - SOra
0 « 3ra/sec
Q • 2.0
oj • 12.0
o, (predicted) • 11.5
x • 100 m
LAGRAN6IAN INTEGRAL
TIM! SCALES:
3 sec
3 sec
0.1
0.2
<*>'
0.3
0.4
0.5
c«v
FIGURE 7. COMPARISON OF ROADWAY MODEL PREDICTED CROSS-WIND
INTEGRATED CONCENTRATION WITH GAUSSIAN PLUME
FORMULA PREDICTIONS AT A DISTANCE OF TOO m
DOWNWIND OF A POINT SOURCE AT 50 m ELEVATION
24
-------
30
ers
t>0
o
10
1
I
• MODEL PREDICTION
* CALCULATED VALUES FROM
GAUSSIAN PLUME FORMULA:
= 5 m
= 3 m/sec
= 2.0
= r. = 3 sec
•X
t
I
1
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.0 0.9 1.0
= c~6v
FIGURE 8. COMPARISON OF ROADWAY MODEL PREDICTED CROSS-WIND INTEGRATED
CONCENTRATIONS WITH GAUSSIAN PLUME FORMULA PREDICTIONS AT
A DISTANCE OF 100 m DOWNWIND OF A POINT SOURCE AT 5 m ELEVATION
25
-------
plume formula predictions is evidence of the validity of our diffusion
modeling approach inasmuch as the plume formula predictions are correct
under the conditions used in this simulation.
26
-------
IV METEOROLOGICAL INPUTS TO THE MODEL
Unlike the GM study in which extensive measurements were made of tur-
bulence quantities, most applied pollution dispersion studies must rely on
a limited set of available observations. For this reason, rather than use
the GM turbulence data directly in the roadway model, we utilize
theoretical-empirical expressions for the quantities required that can be
evaluated given only the local values of surface roughness z0, Obukhov
length L, and mean wind speed u,n at 10 m. This chapter discusses these
expressions in detail and compares the values they yield with those from
the GM wind and turbulence data.
In Chapter II, we indicated that the model requires profiles of the
jy
mean wind component u(z) normal to the roadway, the profile u (z] of the
same component of the turbulent velocjjty_ variance, the profile of the
vertical turbulent velocity variance w'^(z), and the Lagrangian integral
time scale T of the turbulence. The last quantity is difficult to measure
because it cannot be obtained directly from turbulence velocity observa-
tions made at fixed points. Draxler (1976) has inferred approximate values
for T from diffusion data collected at a number of sites. His estimates
of T are of the order of 100 sec for the horizontal velocity component
under both stable and unstable conditions. His suggested values for the
time scale of the vertical velocity fluctuations are 60 sec under unstable
conditions, and 30 sec for stably stratified flows. We found that the
model predictions agreed better with the observations at the most remote
sampling sites (where concentrations are most sensitive to T} if integral
time scale values of only one-half those suggested by Draxler were used
for the vertical velocity fluctuations. Thus, in all of our simulations,
we used
27
-------
= 100 sec
30 sec , if L < 0 ;
lw
10 sec , if L > 0 ,
where L is the Monin-Obukhov length.
— 2 2
To formulate the profiles of u, u' , and w' , we must first express
the direction of the mean wind in the coordinate system shown in Figure 9.
Recall that in our model u and u' are wind components normal to the road-
way, regardless of the orientation of the road. Keep in mind, too, that
the component 7 of the wind parallel to the highway enters only the source
strength term [see Eq. (8) or Eq. (9)] and that the corresponding turbu-
lence component v' does not enter at all. This occurs because of our
assumption of symmetry of the roadway with respect to the receptor plane.
The expression for u(z) is the following:
u(z) = ^g(z,L,zQ)sin e , (22)
where u* is the friction velocity, k = 0.4 is the von Karman constant, L is
the Monin-Obukhov length, ZQ is the surface roughness, and e is the mean
wind direction in the roadway coordinates shown in Figure 9. The form of
the function g is determined from similarity theory and depends on the
sign of L as follows [fromflagland (1973)]:
/z + zn\
1/L = 0: g = £n (— - -^ ; (23)
(24)
\ * - ]N*n + ])
) + £"-M)(*-1) ;
L < 0: g = 2 tan' * - tan" *J + In /.,. . lW... . 1J ; (25)
28
-------
where
4* = 1 - 15
(26a)
*o= M5r
(26b)
ROADWAY-
RECEPTOR PLANE
(x-z) LOCATION
FIGURE 9. COORDINATE SYSTEM USED TO EXPRESS MEAN WIND
AND TURBULENCE VELOCITIES
Th.e expressions for 1/^(2) and w'2(z) are derived from the semi-
empirical formulation of Binkowski (1978). He gives expressions for the
rms vertical velocity fluctuation o and the fluctuating components o
and o normal and parallel, respectively, to the direction of the mean
wind. Translating these formulas into the coordinate system of the road-
way shown in Figure 9, we obtain
where
(27)
u'2(z) = [cos 2F]"1 • [
oj +
cos2e - o +
(28)
(29a)
29
-------
2 2 2^2
o = Q - J • o
U ^ W V
?r
= <£ 3.
L
0 + 0.75
m
1.8r
(29b)
(29c)
rm
- L
jo - i6,r1/4 ,
(l + 5; , ;> 0 ,
0.4[0.4 + 0.6 exp(4?)j
0.4(1.0 + 3.4? - 0.25;2)
o
(29e)
(29f)
c < o
(29f)
z + z,
(29g)
and e is the angle of the mean wind direction shown in Figure 9.
Equations (22) through (26) constitute the bases of the mean wind esti-
mates used in the roadway model, and Eqs. (27) through (29) make up the tur-
bulence velocity variance equations employed. Both sets of equations require
as inputs the friction velocity u*, the Monin-Obukhov length L, and the sur-
face roughness. The remainder of this chapter describes (_1_) our method of
checking the agreement between the values of IT, u'2, and w'2 given by Eqs.
(22) through (29) and those measured during the GM study, and (2) the deri-
vation of an empirical formula from which estimates of u* and L can be
obtained using wind speed data alone.
To obtain u* and L values for each of the GM experiment periods, we used
u^ and L as fitting parameters to arrive at least-squares fits of the mean
wind profiles given by Eqs. (22) through (26) to the GM tower measurements.
In the process, we also obtained a value for the surface roughness of the GM
site. The following procedure was used.
30
-------
First, we chose a value of ZQ in the range of values that is reason-
able for the terrain surrounding the GM test facility, i.e., 0.01 < zn < 0.1
meter. Then, consideration of the wind data from Towers 1, 5, and 6 only
(the other towers are so close to the roadway that they are subject to vehicle
wake influences) resulted in a total of nine observations of u\ three at
each of the levels of 1.5, 4.5, and 10.5 meters, to which the least squares
fitting routine could be applied. Letting e be the difference between the
observed value of u" at station n (n = 1, 2, ..., 9) and the value calculated
from the formula u = u*g(zn,L,zQ)/k, we chose u* such that the mean error,
T, was zero. The following equations were written:
" . U*9n . n = 1, ..., 9 ; (30)
en n " k
where
u = observation at station n,
gn - g(zn.L,z0)
and
n
, = 0
or
for given ZQ and L.
Thus,
(3D
A computer program was set up to solve the u* values for each experiment
from the initial guess of the surface roughness (ZQ). Using the following
equation, ICO Monin-Obukhov lengths were generated:
31
-------
L = exp(X) , (32)
where
X = 0, 0.1, 0.2, ..., 10.0
and
sic]n(L) = sign(Ri)
where R. is the Richardson number measured during each experiment. For each
of the 100 values of L, we determined u* from Eq. (31) and the corresponding
theoretical value ~u at each of the nine observation sites from Eq. (22), namely,
u{z) = ^fg(z,L,z0) . (33)
~~2
The mean square difference e (L) between the observed wind speeds and
those given by Eq. (32) for each of the 100 selected values of L was deter-
mined by:
9
2 2
c(L) - [un - u(zn,z0,L)] . (34)
n~* i
Of the 100 values of L, the best estimate of the true value was taken to
~2 ~~7
be that for which e is a minimum. The minimum values of e* for each of
the GM experiments were then summed to obtain a "global" error for the
chosen value of z~. The entire process described above was repeated using
a different guess for z,.; from the set of global errors obtained, the
true value of the surface roughness ZQ was taken to be that corresponding
to the minimum global error.
Using the technique outlined above, we found that a surface roughness
ZQ = 0.05 m gave the best fit of the observed wind data to the similarity
profiles. The values of u* and L that we derived for this surface rough-
ness are listed in Table 1 for each of the GM experiment periods.
32
-------
TABLE 1. ESTIMATED VALUES OF u* AND L FOR THE
GM EXPERIMENT PERIODS (ZQ = 5 cm) AND
THE ANGLE ¥ OF THE MEAN WIND RELATIVE
TO THE ROAD
Qbukhov Length
GM Experiment
275080
275083
275090
275093
276081
276084
276091
276094
279080
279084
279090
279093
281080
281083
281090
281093
283081
283085
283092
283095
286081
286084
286091
286094
290080
290083
290090
290093
293103
~~ (mf "
-2.7
-2.7
-2.7
-2.7
-270.4
-181.2
-99.4
-90.0
16.4
24.5
-121.5
-109.9
73.7
298.9
148.4
1.8
54.6
-20.1
-2.7
-66.7
110.0
73.7
2.2
-330.3
812.4
181.3
601.9
1096.6
-544.6
(m/sec)
0.20
0.29
0.37
0.41
0.20
0.20
0.27
0.30
0.09
0.09
0.13
0.18
0.17
0.11
0.08
0.19
0.11
0.10
0.15
0.13
0.21
0.23
0.17
0.25
0.19
0.18
0.19
0.22
0.19
(rad)
2.50
2.28
2.76
2.81
0.56
0.64
0.83
0.99
1.23
1.14
1.28
1.23
3.96
3.86
4.76
4.85
1.25
1.39
1.01
0.64
0.19
0.33
0.41
0.45
4.38
4.15
4.20
4.11
1.58
33
-------
TABLE 1 (Concluded)
Obukhov Length
GM Experiment
293110
294080
294083
294090
294093
295080
295083
295090
295093
296080
296083
296090
296093
297080
297083
297090
300080
300083
300090
300093
302080
302083
302090
302093
303080
303083
303090
303093
(m)
-330.3
18.2
21.0
-200.3
-66.7
4.1
18.2
3.7
2.7
134.3
181.3
244.7
-134.3
99.5
90.0
200.3
36.6
-121.5
-200.3
-99.5
1.8
-30.0
-36.6
-30.0
1.5
-1.1
-1.1
-5.0
(m/sec)
0.19
0.14
0.14
0.10
0.13
0.01
0.04
0.03
0.02
0.23
0.25
0.28
0.26
0.20
0.19
0.25
o.in
0.19
0.21
0.21
0,19
0.23
0.28
0.27
0.11
0.11
0.14
0.21
(rad
1.57
0.87
0.86
0.63
0.77
4.65
4.25
4.15
3.70
0.04
0.04
0.16
0.19
0.06
0.07
6.24
0.28
0.43
0.37
0.41
2.84
2.97
3.06
3.15
2.47
2.97
3.34
3.54
34
-------
Table 1 also includes the angle e of the mean wind. Figures 10(a) through
10(c) compare the observed wind speeds at each of the three levels of 1.5,
4.5, and 10.5 meters with the corresponding values obtained from Eqs. (22)
through (26) using the u* and L values shown in Table 1. The observations
represent average wind speeds at each level in Towers 1,5, and 6. A vast
majority of the predicted speeds are within ±10 percent of the measured
values, indicating that the similarity expressions incorporated in Eqs. (22)
through (26) are indeed accurate representations of the wind speed profiles
near the ground.
Checking the accuracies of the turbulence velocity variance expres-
sions [Eqs. (27) through (29)] requires some caution because the measured
values of these parameters are more sensitive to vehicle wake turbulence
than are the mean wind components. Accordingly, we first considered only
situations with westerly winds, i.e., wind directions between 250° and 290°,
and used only the turbulence data collected at Tower 1, which is farthest
west of the roadway. The resulting observed and predicted values of
u' = ~^Z * and w1 = w'2 ' are displayed in Figures ll(a) and ll(b),
respectively. The observed and predicted values of w1 agree quite well--60
percent of the predicted values are within ±20 percent of the measurements;
but the predicted values of u1 are consistently too low by about 35 percent.
Almost all of the data points used in Figure 11 are for cases of unstable
stratification, and it is during these conditions that the Kansas and
Minnesota turbulence data, which Binkowski (1970) used to develop the
formulas in Eq. (29), show the largest scatter in the o and o values.
In fact, the scatter is sufficiently large that profiles of o and a
exist that would be compatible with both the Binkowski and the GM data
sets. We have not pursued this point. In the model validation studies
(discussed in Chapter V), we simply multiplied all the u1 predictions by
1.5 to eliminate the discrepancy shown in Figure 11 (a).
35
-------
i.t
a.» a.*
Observed
4.*
I.t
(a) z = 10.5 meters (zn = 0.05 meters)
FIGURE 10. OBSERVED VERSUS PREDICTED MEAN WIND SPEED
36
-------
3.0 -
Cl
4->
U
•n
c
i .0
3 »
Observed
0.*
(b) z = 4.5 meters (z = 0.05 meters)
FIGURE 10 (Continued)
37
-------
o
4->
'_>
3.*
Observed
(c) z = 1.5 meters (ZQ = 0.05 meters)
FIGURE 10 (Concluded)
38
-------
0.0
0.2
0.8
-0.8
r.o
0.2
0.4 0.6
Observations
0.8
0.0
—pl/2
(a) u'2
•1/2
FIGURE 11. OBSERVED VERSUS PREDICTED VALUED OF u'2 AND w'2
Only cases with wind direction e in the range
250 < ¥ < 290 are represented, and only data from
Tower 1 are included.
1/2
39
-------
0.0 0.1
±20% ERROR BOUNDS -
0.2 0.3
Observations
(b) w
,2
1/2
FIGURE 11 (Concluded)
40
-------
Our efforts to check the accuracies of the variance formulas [Eqs. (27)
through (29)] were extended by considering data from Tower 6, the east-most
tower, only in cases of wind directions e" in the range of 70° < ¥ < 110°.
The predicted and observed values of u1 and w1 are shown in Figures 12(a)
and 12(b), respectively. The tendency for the predicted u1 value to be too
low is again apparent. One noteworthy set of points shown in both figures
is that for which the predicted values of u1 and w' are near zero--u' = 0.02,
w' = 0.01--but the observed values are over 10 times larger. This particular
set of points corresponds to a case of extremely stable flow and light winds
(Day 295) with a vertical temperature gradient of +1°C per 5 meters and
u ^0.4 m/sec. The large value of u1 measured in this case is most likely
due to meanders in the local circulation induced by topography and inhomogen-
eities in surface roughness, and the large w' values observed may be due to horv
zontal divergence in the circulation, to gravity waves, or to both. In any
event, the fluctuating velocities observed are not the result of the mech-
anical and buoyancy phenomena that generate the turbulence described by
similarity theory. Describing the dispersive properties of very stably
stratified flows is one of the largest unsolved problems of diffusion
modeling.
Finally, we compared the predicted u' and w1 values with those mea-
sured at Towers 1, 5, and 6 when the wind direction e" fell in either of
the sectors 340 < ? < 20 or 160 < ¥ < 200. These comparisons are given in
Figures 13(a) and 13(b). The dashed line in Figure 13(a) is the best fit
line of the u' values shown in Figure ll(a). It is apparent that with
northerly or southerly winds the predicted u' values are not as low as
they are for westerly flows. The dependence on wind direction of the dis-
crepancy between the observed and predicted values of u' may be the result
of the dependence of o and a on inversion height. Panofsky et al. (1977)
have proposed the relationship:
ou = ay = u*[12 - h/(2L)]1/3 , L < 0 , (35)
41
-------
0.0
0.2
0.6
0.8
0.8
0.8
0
•5
cu
Q.
O.B
0.4
O.B
0.2
0.2
0.2
J I
.4 0.6
Observations
0.8
0.0
FIGURE 12.
(a) u' = u'2
1/2
OBSERVED VERSUS PREDICTED VALUES OF u' = u'2
1/2
AND w' = w'2
1/2
Only cases with wind direction e in the range 70 < e < 110 are
represented, and only data from Tower 6 are included.
42
-------
0.0
0. t
0.1
0.2 0.3
observations
1/2
-0.1
D.O
(b) w1 = w'2
FIGURE 12 (Concluded)
43
-------
0.0
0.2
0.4
0.6
0.8
•*->
o
•^
-o
OJ
Q_
O.J
0.2
0.4 0.6
Observations
(a) u'
0.8
0.0
FIGURE 13. OBSERVED VERSUS PREDICTED VALUES OF u1 AND w1. Only
data for cases with wind direction in the sectors
34Q < Q < 20 and 160 < "5" < 200 are represented, and
only measurements from Towers 1, 5, and 6 are included.
44
-------
0.0
0.1
0.2
0.3
o.q
O.J
0.1
0.2 0.3
Observations
(b) W
0.1
0.0
FIGURE 13 (Concluded)
45
-------
where h is the mixed layer depth or inversion base height. They have found
very good agreement between this equation and measured data. Binkowski's
formulas [Eq. (29)] for o and a are independent of h. Consequently, it
is quite likely that the dependence on wind direction of the accuracy of
Eq. (29) is a result of a systematic variation of h with wind direction at
the GM site. Since inversion data were not measured during the GM study,
we cannot test this hypothesis now.
On the whole, the comparisons presented in Figures 11 through 13 indi-
cate that Eq. (27) provides acceptably accurate estimates of w'2 ' but
_ -I /O
that the u'2 ' values obtained through Eq. (29) are subject to some uncer-
tainty. We believe that through proper modifications of Eq. (29) or by
replacing this formula altogether with Eq. (35), acceptably accurate esti-
mates of the horizontal turbulence variances can also be obtained.
As a result of the good agreement between the if values predicted by
similarity theory and the measurements (see Figure 10), we decided to
explore the possibility of deriving estimates of u* and L from wind data
alone. According to Eq. (22), the ratio of the wind speeds measured at
two heights, z, and Zp, is:
IT, g(z,,L,zn)
_ l_ _ I __ U__
TI2 " g(z2,L,zQ) ~ '
For given values of z, , z^, and ZQ, the right side of Eq. (36) is a
unique function of L, at least in the interval -0.1 <. 1/L <. 0.1. Figure 14
shows a plot of this function for z, =" 10.5, z? = 1.5, and z~ = 5 cm. Also
shown is a plot of the drag coefficient:
c = -u-* =
CD
obtained from Eqs. (22) through (26). Using the wind speed measure-
ments made at 10.5 and 1.5 meters and the two curves plotted in Figure 14,
46
-------
LT)
O
I =3
*
0.10 -
0.08-
0.06 -
0.04 -
-0.10
-0.05
0.05
1/L (per meter)
= 2.5
-2.0
LO
O
-1.5
1.0
0.10
FIGURE 14. PLOTS OF THE RATIO OF THE WIND SPEEDS AT 10.5 AND 1.5 METERS
AND THE DRAG COEFFICIENT u*/u(10.5) GIVEN BY SIMILARITY
THEORY AS FUNCTIONS OF 1/L. The surface roughness, ZQ, is 5 cm.
-------
we can obtain estimates of both u* and L. For convenience, we refer
to this as the wind speed ratio (WSPx) method. The technique can be used
for any values of z,, z2, and ZQ by recomputing the Zu/ZL and u*/u., curves.
We should also point out that with this method, L = 10 is assumed if Lu/IL
exceeds the theoretical value at 1/L = 0.1, and L = -10 is assumed if the
wind speed ratio is less than the theoretical value at 1/L = -0.1. Figure
15 compares the u* and L values given by the WSR method for each of the GM
experiments with the corresponding values listed in Table 1. In applying
the WSR method, we used data from only Tower 1 or Tower 6, depending on
which was upwind of the roadway, rather than averages for several towers.
Looking first at the u* estimates, we find that the two sets of values
are in fair agreement: Over 67 percent of the WSR estimates are within ±20
percent of the values given in Table 1. However, there is a distinct ten-
dency for the WSR method to overpredict u*, a direct result of the method's
tendency to predict unstable stratification {for which u*/lT,) is large) in
situations where the values in Table 1 indicate that stable stratification
existed. Recall that in obtaining the L values of Table 1, the sign of L
was taken from the Richardson number, which was measured during each of the
GM experiments. In general, there is little agreement between the two sets
of 1/L values shown in Figure 15(b). In fact, only 60 percent of the values
given by the WSR method agree in sign with those in Table I—only
slightly better than the agreement that could be obtained by making random
guesses of the sign of L.
A case-by-case examination of the meteorological conditions existing
during those periods when the signs of the L values given by the two
methods disagree revealed the following findings. On 4 of the 23 occa-
sions for which the signs differ, the sign of the Richardson number mea-
sured at Tower 1 was opposite that recorded at Tower 6 some 100 meters
away. In 9 of the remaining 19 cases, the temperature stratification was
not uniform vertically at the tower from which the wind data used in the
WSR method were taken; the sign of the temperature difference between the
top two levels of the tower was opposite that between the bottom two levels.
The 4 cases in which the Richardson numbers at Towers 1 and 6 had opposite
48
-------
+15%
0.3
o
.c
-i->
-------
-0.1 -0.08 -0".06 -0'.04 -0'.02'"V
.* "0.02 0.04 0.06 0.08 OJO
Observed 1/L (Table 1)
-0.02
• 4
-0.04 K
3
V-0.06 -
4 -0.08
+ -0.1
•o
0)
4->
o
•I—
•o
o
CL
(b) Monin-Obukhov Length, L
FIGURE 15 (Concluded)
50
-------
signs all occurred on Day 290, when the flow was stably stratified upwind
of the roadway during the observation period and unstable downwind. This
suggests that the heat from large numbers of automobiles and the release
of heat stored in the roadway itself can significantly alter the turbulence
locally when the ambient flow is stably stratified. Another possible cause
of change in stability downwind of the road is the overturning of the cold
air near the ground and the warmer air above by the vehicle wake turbulence.
Our finding that the vertical temperature structure near the ground was
not uniform in a number of instances in which the two methods of estimating
L give radically different results emphasizes the problem of characteriz-
ing stability when the vertical temperature profile is nearly isothermal.
Chapter V quantitatively examines the magnitude of the errors in mean con-
centration estimates resulting from uncertainties in the values of u* and L,
51
-------
V RESULTS OF PRELIMINARY MODEL
VALIDATION EXPERIMENTS
This chapter discusses the SFg concentrations predicted by our model
for each of the GM experiment periods and compares them with the values
actually observed. The roadway used in the GM study is oriented north-south.
Instrumented towers were placed 30 and 2 meters west of the roadway; one
tower was sited in the roadway median; and additional towers were located
at distances of 4, 15, 30, 50, and 100 meters east of the roadway's eastern
edge. Thirty-minute averaged SFg concentrations were measured at the 1.5,
4.5, and 10.5 meter levels of all towers except those located 50 and 100
meters east of the road, where measurements were made at 1.5 meters only. Thus,
a total of 20 SFg measurements were made during each experiment. The reader
is referred to Cadle et al. (1976) for further details on the GM monitoring
network and study conditions.
A. SIMULATIONS THAT NEGLECT WAKE EFFECTS
In previous chapters we detailed all aspects of the roadway model
except one—the vehicle wake turbulence parameterization. Most early
roadway models ignored the dispersive effects of wake turbulence entirely,
but models developed within the last few years have included some means
of approximating these effects. Apparently, however, no quantitative model-
ing studies have ever been performed to establish the actual importance of
wake effects to accurate predictions. We performed such a study by
exercising our model under the conditions of the GM experiments, with a
wake turbulence parameterization omitted entirely. Obtaining good agree-
ment between the simulation results and the measured data would indicate
that, at least insofar as the mean concentration of an inert material (SFg)
is concerned, wake effects are minor and can be neglected.
52
-------
In Figure 16, the SFg concentrations predicted by the model with a
wake parameterization omitted are compared with the GM data. We emphasize
that the only input data used in the model for these calculations were the
SFg emissions rates reported in Cadle et al. (1976), the u*, F, and L values
given earlier in Table 1, and estimates of the time scales T and T given
on page 28. The only exceptions are the four experiments conducted on Day
295. As we noted earlier in Chapter IV, similarity theory completely failed
to duplicate conditions on this extremely calm, stable day. The bulk
Richardson's numbers measured during the four experiments conducted on
Day 295 ranged between 0.35 and 7.5, which are outside the limit (viz, 0.2)
of applicability of the similarity expressions that we use in the model.
Consequently, in the simulations of the four experiments conducted on
Day 295, we forced the IT, a , and o profiles to agree with the averaged
measured values at the 10.5 meter levels of the towers 30 meters east and
west of the road. Our justification for altering the input data in this
manner lies in the fact that we are concerned here with how well the dif-
fusion model itself performs under various meteorological conditions.
Therefore, we are attempting in these validation and testing exercises to
minimize errors in the input data as much as possible so that errors in
the predicted concentrations will reveal inadequacies in the model itself.
We showed in the previous chapter that the similarity expressions that we
have incorporated in the diffusion model yield quite accurate mean wind
and turbulence profiles given only the local values of the friction velocity
and Monin-Obukhov length. However, under extremely stable conditions,
namely R. > 0.2, these expressions break down altogether and the profiles
they produce are completely unreliable. Consequently, under these con-
ditions one must rely on measured winds and turbulence parameters for
reliable estimates. We have adhered to this procedure in our validation
studies.
We find in Figure 16 that, with the exception of several of the values
predicted for Day 295 (designated by arrows in the figure), the performance
of the model with wake effects ignored is not particularly bad. The rms
error of the 1100 predicted values is 0.80 ppb (the averaged observed value
53
-------
12.5
en
ii.i
7.5
on
u
-------
is 0.71), the correlation coefficient p = 0.79, and the regression line is
PRED = 0.07 +1.27 • OBS. One of the features that sets Day 295 apart from
all the others is extremely light winds: IT--0.5 m/sec. On all other days,
wind speeds were typically higher than 1.5 m/sec.
Considering the results shown in Figure 16 from this perspective, we
speculate that wake effects are apparently most pronounced when wind speeds
are less than 1 m/sec, at least when the vehicle speeds are 22 m/sec and
higher, as they were in the GM study. This speculation is supported by
Figure 17, where we have replotted the results shown in Figure 16 with
the 80 values predicted for Day 295 deleted. In this subsample of 1020
points, and with the wake parameterization scheme still omitted, the rms
error is o = 0.53 ppb, the correlation coefficient p = 0.87, and the regres-
sion line is PRED = 0.07 +1.16 • OBS.
B. SIMULATIONS THAT INCLUDE A SIMPLE PARAMETERIZATION OF WAKE EFFECTS
At least three quantities are required to model the wake properly:
the initial turbulent energy in the wake, the time scale of the energy
decay, and the Lagrangian integral time scale of the wake turbulence.
Because values are not currently available for any of these parameters,
it was necessary to formulate a wake parameterization empirically.
Our initial efforts, which are quite primitive, were based on the follow-
ing assumptions. First, for simplicity, the integral time scale of the
wake turbulence was assumed to be zero. (This is tantamount to assuming
that the diffusive action of the wake is akin to that of Fickian diffusion.)
Second, the initial wake turbulence energy was assumed to be proportional
to the work done by the drag force on the motor vehicle; that is,
55
-------
it
en
CTl
-Q
§:
oo
-a
0.25) OMITTED
-------
(38)
where
Cn = the vehicle drag coefficient,
A = the cross-sectional area of the vehicle,
V = the vehicle speed relative to the air,
y ,y = constants,
2 2
u" ,w" = the horizontal and vertical components,
respectively, of the initial mean square
velocity fluctuations in the wake.
Third, to account for the decay of the wake energy, the wake intensity
was assumed to be constant at the level given by Eqs. (38) and (39) for a
period of T seconds, after which the energy was assumed to drop abruptly
to zero and to remain there for all later travel time.
The foregoing description of the wake was incorporated into the
stochastic velocity generator used to produce the simulated wake velocity
components u" and w" described in Chapter II [see Eq. (15)].
To avoid using y , y , and T as devices for "tuning" the model, we
U W W
attempted to determine values for the first two parameters from the wind
and turbulence measurements made at the sites closest to the roadway.
We chose
•y = Y
u rw
2
and then using CD = 0.5, and A = (2.5 m) , values chosen somewhat arbi-
trarily, we selected a value for y that made Eqs. (38) and (39) yield
velocity variances close to those suggested by the wind and turbulence
57
-------
measurements. Through this process, we found 0.005 < y < 0.01, Next,
to refine the estimate of y and to obtain a value for T , we exercised
the model with the wake parameterization included to simulate Experiment
279080. (This period was chosen completely arbitrarily.) We varied the
value of y within the range noted above, and we experimented with various
values of T to determine which set of values produced concentration
w
predictions in closest agreement with the measured data. We were partic-
ularly concerned with reproducing the upstream diffusion that did not
occur in the simulation in which the wake was ignored. The following
parameter values were selected from these exercises:
Tw - 3 sec ,
Y = 0.007 sec1/2/m3/2
Using these values in the wake parameterization along with the values
of C~ and A given earlier, we incorporated the wake simulation algorithm
into the model and then repeated the simulation of all of the GM experi-
ment periods. The source emission rates and the values of u*, 9", L, T , and
T employed in these calculations were identical to those used in the simu-
lations illustrated in Figure 16 in which the wake effects were ignored.
Figure 18 shows the results of the comparison of the new calculations {with
the wake simulation included) with the observed data.
Comparison of the results presented in Figure 16 and 18 indicates that
the wake parameterization scheme, as simple as it is, has improved the
accuracy of the model predictions considerably. The correlation coefficient
of the new results is very high (p = 0.91), and 53 percent of the 1100 pre-
dicted concentrations are within ±30 percent of the measured values (as
illustrated later in Figure 22). The rms error a = 0.33 and the regres-
sion line is PRED = 0.13 + 0.91 • OBS. In the course of testing the wake
parameterization scheme, we discovered that while the overall performance
of the model was greatly improved with the scheme included, there were a
number of simulation periods in which the model performed considerably
58
-------
en
10
2 3
Observed SFg (ppb)
FIGURE 18. COMPARISON OF OBSERVED SF6 CONCENTRATIONS WITH THOSE PREDICTED BY THE
ROADWAY MODEL WITH WAKE TURBULENCE EFFECTS PARAMETERIZED: ALL
EXPERIMENTS (1100 DATA POINTS) INCLUDED
-------
better with the wake parameterization omitted. On seeking an explanation
of this finding, we discovered that in all cases in which the model per-
formed better with the wake effects ignored, the mean wind direction e"
was within about ±10° of the roadway, namely, (sin ¥ | < 0.2. The number
of cases involved (13 of the 55 GM study periods), the high correlation
with wind direction, and the difference in the rms error of the model
realized between simulations made with and without the wake parameteriza-
tion were strong evidence that the phenomenon in question was not merely
a chance event, but rather the result of some physical process associated
with the vehicle wake.
After examining the condition of the GM experiment more closely, we
concluded, with some confidence, that the process in question is the so-
called "slipstream" phenomenon. The vehicles in the GM experiment were
spaced only a few vehicle lengths apart, so that each vehicle was within
the wake of the one ahead of it, except when the mean wind component normal
to the roadway was strong enough to force each vehicle's wake out of the
path of the vehicle immediately behind. When the vehicles are well within
each other's wakes, the air speed relative to each vehicle is considerably
lower than that which a single vehicle would experience in the same meteo-
rological conditions. Thus, with close vehicle spacing, the drag work done
by each vehicle, and hence the intensity of turbulence in the wake pro-
duced, are lower than that which would occur if each vehicle were not
within the wakes of others.
To account for this phenomenon, we set the air velocity relative to
each vehicle to 1/10 of its nominal value when the wind angle e" is in the
range |sin ¥| < 0.2. The model results displayed in Figure 18 include
this modification.
In order to analyze the model's capabilities further, we have plotted
in Figures 19 and 20 scatter diagrams of predicted versus observed concen-
trations for stable and unstable conditions separately. These figures
reveal a tendency for the model to overpredict slightly in stably
60
-------
Observed SFg (ppb)
FIGURE 19. COMPARISON OF OBSERVED SF6 CONCENTRATIONS WITH THOSE PREDICTED BY THE
ROADWAY MODEL WITH WAKE TURBULENCE EFFECTS PARAMETERIZED: RESULTS
DISPLAYED FOR STABLY STRATIFIED FLOWS ONLY (580 DATA POINTS)
-------
ro
I I I | I I I I 1 I I I 1
l.S l.t 8.5
Observed SF6 (ppb)
3.S
FIGURE 20. COMPARISON OF OBSERVED SF6 CONCENTRATIONS WITH THOSE PREDICTED BY THE
ROADWAY MODEL WITH WAKE TURBULENCE EFFECTS PARAMETERIZED: RESULTS
DISPLAYED FOR UNSTABLY STRATIFIED FLOWS ONLY (520 DATA POINTS)
-------
stratified flows and to underpredict during unstable situations. These
tendencies are also apparent from the following performance statistics:
> Stable cases only
Number of data points N = 580
Average predicted concentrations P = 0.85 ppb
Average observed concentrations 0 = 0.76 ppb
a = 0.38 ppb
P = 0.91
PRED = 0.13 + 0.95 • OBS.
> Unstable cases only
N = 520
P = 0.68 ppb
0 = 0.65 ppb
a = 0.28 ppb
p = 0.91
PRED = 0.14 + 0.83 • OBS.
The fact that the intercept of the regression line for the cases
with the wake scheme employed is larger than that obtained with the scheme
omitted indicates that the wake parameterization is spreading material
over too large an area. Another feature of the model results is the lower
slope of the regression line achieved. Contributing to this deficiency is
the set of points predicted for Day 275. For some unknown reason, the
highest observed values on one of the experiments of this day are consis-
tently higher than the predicted values by a factor of two. The other
experiment shows a similar underestimation.
A quantity of great concern in applied air pollution studies is the
expected peak concentration. Figure 21 displays a scatter plot of the
predicted versus observed peak values obtained in the 55 experiments. The
sloping dashed lines in the figure delineate the ±25 percent error bounds.
Over 75 percent of the calculated peak values fall within these bounds.
63
-------
PRED = -
1.25 • OBS
CL
10
u_
CO
"So 3
o
• r-
•o
(!)
Q-
RED =
0.75 • OB
III! 'I 1 I I I I I
1 1 1 »
fill
1.6
2.*
8.6 3.« 3.6 «.»
Observed Peak SFg (ppb)
4.S
6.»
FIGURE 21. COMPARISON OF OBSERVED PEAK SFe CONCENTRATION FOR EACH OF THE
55 GM EXPERIMENTS WITH THE VALUES PREDICTED BY THE MODEL
(WAKE EFFECTS INCLUDED)
-------
The performance statistics of the model are summarized in Table 2,
and plots of the fractional error distributions found for several of the
validation studies are given in Figure 22. The basic model requires about
36,000 words of core storage and it used about 10 minutes of CPU time on
the UNI VAC 1110 to simulate all 55 of the GM experiments. Shorter execu-
tion times can be achieved by relaxing the error bounds on the predicted
concentrations (see Table 2). The computer memory requirement increases
by 3 words for each additional monitor site added, regardless of the dis-
tance of the monitor from the roadway. The machine time required by the
model increases very slightly as additional sampling sites are added to
a given area. However, if the area itself is extended, the machine time
requirement increases. We have not established the quantitative dependence
of the machine time requirement on the size of the sampling domain.
65
-------
TABLE 2. MODEL PERFORMANCE STATISTICS
Error Distribution--
percentage of
predicted value
with error <.:
Regression Line--
PRED = a + B -DBS
1.
2.
3.
4.
5.
N =
P" =
0 =
a =
P
Conditions
All experiments,
wake "on"
Day 295 omitted,
wake "on"
Same as No. 1 but
with relaxed
error (execution
time - 3 min;
Case 1 = 10 min)
All experiments,
wake "off"
Day 295 omitted,
wake "off"
total data points.
Average predicted SF
Average observed SFg
rms error.
N P 0 o p 25% 50% 100% a
1100 0.76 0.71 0.33 0.91 47% 68% 83% 0.12
1020 0.74 0.70 0.30 0.92 48 70 85 0.11
1100 0.76 0.71 0.38 0.88 42 65 81 0.13
1100 0.96 0.71 0.80 0.79 36 57 80 0.07
1020 0.88 0.70 0.53 0.87 37 59 82 0.07
6 (Ppb).
•
e
0.91
0.90
0.90
1.27
1.16
correlation coefficient.
-------
100
cu
C1 C
(D •!-
-l-> -C
C 4J
0> •..-
U 3
i-
OJ
ex
80
CO
cu
O
Si 60
u
•»- C
-a cu
a> >
40
20
WAKE "ON"
WAKE "OFF"--DAY 295
EXCLUDED
WAKE "OFF"
20
40 60
Error (percent)
80
100
120
FIGURE 22. DISTRIBUTION OF FRACTIONAL ERROR. The curve labeled "wake 'on'" includes
wake parameterization; "wake 'off'" simulations ignore wake effects.
-------
VI SUMMARY
We have developed a model of air pollutant transport, dispersion,
and removal processes based on Lagrangian principles. The model requires
as inputs estimates of the local friction velocity, the Monin-Obukhov
length, the surface roughness, and source emission rates. In applications
to pollutant dispersal within 100 meters of a roadway, we showed that the
model performs very well provided that the Richardson's number does not
exceed approximately 0.25. We also found that
> Vehicle wake turbulence plays a dominant role in
dispersing material within 10 meters of roadways when
the wind speeds are less than about 1 m/sec and the
vehicle speeds are about 20 m/sec or larger.
> The dispersive effects of wake turbulence are
greatly reduced when vehicles are closely spaced and
the wind is parallel to the roadway.
> Except for situations in which wind speeds are less
than 1 m/sec, spatial and temporal variations in
material concentration are dominated by ambient tur-
bulence. A correlation coefficient of 0.87 was obtained
under these conditions with wake effects ignored
entirely. However, the rms error of the estimates is
greatly improved if wake effects are taken into account.
68
-------
APPENDIX A
A SCHEME FOR SIMULATING PARTICLE
PAIR MOTIONS IN TURBULENT FLUID
Robert G. Lamb*
Meteorology and Assessment Division
Environmental Sciences Research Laboratory
Environmental Protection Agency
Research Triangle Park, North Carolina 27711
Preprint. Submitted for publication to the
Journal of Computational Physics, 1978
* On assignment with the EPA from NOAA, U.S. Department of Commerce.
69
-------
APPENDIX A
A SCHEME FOR SIMULATING PARTICLE
PAIR MOTIONS IN TURBULENT FLUID
ABSTRACT
A technique is developed and tested for simulating the random
motion of particles in turbulent fluid. The technique is applicable to
single particles or to pairs of particles whose motions are correlated.
The method is intended for use in determining one- or two-particle
Lagrangian turbulence statistics from instantaneous velocity fields
generated by numerical turbulence models. The technique can also be
used as the basis of the Monte Carlo type of models of turbulent diffu-
sion. It requires as inputs the mean, variance, and Lagrangian inte-
gral time scale of the turbulence, and for particle pair simulations,
it additionally requires the energy dissipation rate of the turbulence.
1. INTRODUCTION
In numerical studies of material dispersion in turbulent fluid,
only a portion of the spectrum of fluid motions can be described in a
deterministic manner. These motions comprise the largest scales of
fluid motion, and their combined effect is represented by the so-called
"mean" velocity in the calculation process. The remaining, small scales
of motion are regarded as "turbulence" and are treated in various ways.
The most common of these in air pollution modeling studies is the eddy
diffusivity hypothesis. In this approach, the effects of small-scale
velocity fluctuations on the mean material concentration c are param-
eterized by a semi-empirical eddy diffusivity K in an equation of
the form:
70
-------
|| + y-vc = V-Kvc + S , (A-l)
where v is the "mean" wind vector referred to above and S is the source
strength function.
There are numerous situations in which either accurate expressions
for K are unavailable or the eddy diffusivity hypothesis itself is
inconsistent with the observed character of material diffusion. An
example of the former is the modeling of air pollution on scales of
1000 km or larger. In this case, the "turbulence" comprises eddy sizes
up to about 100 km, i.e., roughly the average separation of wind monitor-
ing stations. An eddy diffusivity representative of this portion of the
velocity spectrum is not available. An example of a situation in which
the eddy diffusivity hypothesis is inconsistent with the character of the
diffusion process is unstable boundary layers. Deardorff and Willis
(1975) have shown that in free convection the vertical component K of K,
which intrinsically is a positive quantity, would have to assume negative
values to describe the diffusion phenomena observed.
The problems of K-theory can be circumvented by working with
Lagrangian diffusion theory. The Lagrangian counterpart of Eq. (A-l)
for c at a specified point x and time t is expressed as follows:
c(x,t) = /YpU.tlxJt'Wxit1) dx1 dt1 . (A-2)
The difficulty with the Lagrangian approach is determining the
form of the kernel p(x,t|xjt')» which is the probability density that a
material particle released at x1 at time t1 will be found at x at time t.
Until recently, the only method of implementing Eq. (A-2) was to assume
a form for p. The Gaussian puff and plume formulas are examples of
models derived from Eq. (A-2) under the assumption that p has the Gaussian
form (see Lamb and Seinfeld, 1973).
71
-------
Within the last few years, numerical turbulence models have been
developed [e.g., Deardorff (1974)] that make it possible to derive
Lagrangian turbulence properties such as p computationally. In these
models, which explicitly resolve large portions of the spectrum of turbu-
lent motions, particles can be released and tracked at will to obtain
the ensembles of particle trajectories required to derive Lagrangian
infornation. Each realization of the ensemble of trajectories needed
to compute p is given by:
rN i
x(t;x0) = I uLxtt'^Qht1] + y'[x(t';x0),t'] dt' + XQ
JQ ' (A-3)
where xn denotes the particle release point (at t = 0); u(x,t) is the
~u ~ ~
resolvable portion of the Eulerian velocity at (x,t), given by the
turbulence model; and u'(x,t) is a random velocity that represents the
combined effect of all fluctuations in the instantaneous velocity field
that the turbulence model cannot explicitly resolve. In general, u'
represents all velocity variations smaller than the model's grid dimen-
sions. This appendix describes a technique for generating stochastic
velocity fields u' suitable for use in calculating Lagrangian turbulence
properties from numerical turbulence data. [For convenience, we shall
hereafter refer to u1 as the subgrid-scale (SGS) velocity.]
Approximation of the kernel p can also be obtained using so-called
Monte Carlo techniques. With these methods (see Chapter II of this report),
p is estimated from an ensemble of particle trajectories created from
computer-generated random velocity fields whose statistics are made to
satisfy given Eulerian statistical properties of the flow, such as velocity
variance profiles, energy dissipation rate profiles, and so on.
For generality, we developed a scheme for generating pairs of cor-
related velocity functions u1,, u'2- Such functions are required to
simulate the trajectories x-,(t;x.|Q) and Xottjx^n) of a pair of particles
72
-------
simultaneously released at points x,g, X^Q located close enough together
that the motions of the particles are correlated for a time. Particle
pair trajectory simulations are needed to determine the two particle
displacement probability density function P2(x,t;x;t'|x,Q,t,;x2Q,t2)
that a pair of particles released at the space-time points (x,n, t,),
- I U I
(XPQ, tp) will be found at (x,t), (xjt')> respectively. The function
p,, is required to calculate the mean square concentration c2(x,t), viz,
t t
c2(x,t) =11 I I P2(x,t;x,t|x;t1;xl,'t"}S(x;t')S(x',1t") df dt" dx' dx"
JJJ0J0 " (A-4)
Shu, Lamb, and Seinfeld (1978) have used this equation in conjunc-
tion with p2 estimates, derived in the manner outlined above from Deardorff's
(1974) numerical simulation of the convective planetary boundary layer, to
model nonlinear chemical reactions among the constituents of a point
source plume and species in the ambient air. Lamb (1976) proposed the
concept of a macroturbulence through which historical meteorological data
can be used to determine a kernel PO, which, when combined with Eq. (A-4)
yields estimates of air pollution episode probabilities.
Deardorff and Peskin (1970), Riley and Patterson (1974), Lamb, Chen,
and Seinfeld (1975), and Lamb (19785) have computed single-particle dis-
placement statistics using Eq. (A-3) and velocity data u from numerical
turbulence models. In all of these studies, the effect of the subgrid-
scale velocity u1 were found to be negligible, and they were generally
ignored. While it is true that the effects of u' on the displacement
statistics of a single particle are small, if u1 constitutes only a
small fraction of the total energy spectrum of fluid motion, the subgrid-
scale velocity component u1 plays a dominant role in determining
the joint displacement statistics of a pair of particles initially
separated by a distance smaller than the dimensions of the grid on which
the model data u are available. This is true even when u1 represents
only a small fraction of the total fluid kinetic energy. For this reason,
the parameterization of the subgrid velocity field is essential in the
73
-------
calculation of mean square concentration distributions. Thus, the main
uses that we foresee for the two-particle scheme described here are in
modeling the mean square concentration by Lagrangian methods, i.e. by
using Eq. (A-4) in conjunction with Eulerian velocity data. The scheme
will also be of value in implementing Eq. (A-2) when the available velocity
data are on such a coarse grid that the subgrid-scale portion of the
energy spectrum is significant; and it can be used as well in conjunction
with Eq. (A-2) as the basis of Monte Carlo diffusion models.
2. DESIGN OF THE PARAMETERIZATION SCHEME
The basic approach to developing an approximation for the subgrid-
scale particle velocities u, and u' is to formulate a rule for generating
random functions such that the statistical properties of those functions
are equivalent to those that the SGS velocities are known to possess.
7
For example, the mean square subgrid-scale velocity u' is usually known.
Thus, the random functions chosen to represent u\ and u~ must satisfy
-^ —j — *- ~i -t.
= = 1
n = u,-, = u1 . The Lagrangian integral time scale of the SGS velocities
and certain mean values of particle displacements that they induce may
also be known. In this case, the selected random functions u-j and ui
must have statistical properties that are consistent with all of these
known quantities. We outline below the SGS velocity and displacement
statistics that are generally known and that serve as constraints on the
formulation of the parameterization scheme. For simplicity, we treat
only the x-component of the SGS velocity and displacement and ignore cross
correlations among the three vector components. Extensions of the scheme
to more general situations are straightforward.
a. Constraints on the Subgrid-Scale Particle Displacements
By definition the mean subgrid-scale velocity u' is zero locally;
and although it does not necessarily follow from this definition that
the mean SGS particle displacements are zero, nevertheless, we assume
that
74
-------
= x]0 , x2(t) = X2Q , (A-5)
where x, Q and X^Q denote the release points of particles 1 and 2.
This assumption will generate an error roughly proportional to the mag-
nitude of the shear in the deterministic flow speed component u at the
points X,Q, X~Q of release. The mean square SGS particle displacement
and separations can be expressed in the following forms:
*(t) = o2(t) + x2Q , (A-6a)
= o2(t) + x2 , (A-6b)
(t)]T= *2(t) + (x1Q - x2Q)2 . (A-7)
- x2
Implicit in Eq. (A-6) is the assumption that x,n and xon are sufficiently
p n ij I U f.\J
close together that o (t) = a0(t) = a (t) . For the moment, we assume
22
that a (t) and s. (t) are known functions.
To specify the random displacements x-, and x2 uniquely, we would
have to specify the entire infinite hierarchy of joint moments of
these variables. In practical situations, at best only the first- and
second-order moments are known. Thus, in developing the SGS parameter-
ization, it suffices to incorporate only Eqs. (A-5) through (A-7)
as constraints.
b. General Form of the Displacement Parameterization
Let x be the displacement of the n-th particle of a pair (n = 1, 2)
produced by our simulation of the turbulence velocity u-' . The constraints
of Eqs. (A-5) through (A-7) on the statistics of x can be satisfied
if xn is of the following form:
75
-------
xn - xnQ = as + bxn > n = 1, 2 ; (A-8)
where a and b are functions of travel time and £ and x are computer-
generated random numbers. To satisfy Eq. (A-5), we must have
xn = 0 . (A-9)
It is also convenient to have 5 and xn statistically independent and Xi
and X? independent of each other. In this case,
cxn = exn = o = 7jx^ . (A-io)
To satisfy Eq. (A-6), we need
and to satisfy Eq. (A-7), we must have
2 - x - x2 = b22 + b22 = a2
(x1 - x2) - (x1Q - x2Q) = bx + bx = a
If we make
= o
2 (A-13a)
a2 , (A-135)
then it follows that Eqs. (A-5) through (A-7) will be satisfied by
Eq. (A-8) if, in addition to Eqs. (A-9), (A-10), and (A-13) we select
a and b such that
76
-------
a2 = i -h?/2o2 I , (A-14a)
b2 = 1 - a2 . (A-14b)
c. General Form of the Velocity Parameterization
The algorithm for generating u' is obtained by differentiating
Eq. (A-7):
u; = a| + ca + bxn + X{f> . (A-15)
With a and b given by Eqs. (A-14a) and (A-14b), respectively, and £; and
x satisfying conditions in Eqs. (A-9), (A-10), and (A-13), Eq. (A-15)
constitutes the basis of our parameterization of the subgrid-scale
2 2
velocity field. Next, we must specify the forms of £, x > o , and £ .
In view of conditions in Eq. (A-13), our choice of forms for c and
x will be guided by our knowledge of a, which until now we have regarded
as a given function. In actuality, we have very little information about
the form of o. We know only that for small times,
a2(t) = u'2t2 , small t , (A-16)
—2
where u1 is the mean square subgrid-scale velocity fluctuation and t is
travel time. The behavior of a at large times is not known. Lamb (1971)
showed that if the turbulence is stationary and the Lagrangian velocity
spectrum is a band-pass spectrum, similar to that which the Eulerian
spectrum of the SGS velocity is known to be, then in the limit of large
travel times the mean square particle displacement a approaches a constant
value. This is in marked contrast to the effect that a white noise
2
velocity spectrum would produce, namely cr ~. t in the limit of large travel
times. Whether the SGS velocities produce a constant, a linearly increas-
2
ing, or an altogether different dispersion o in the limit of large travel
77
-------
times is not known, but we expect that the details of the variations in
the long-time regime are not crucial. Even if particles are released near
a boundary where the SGS velocities constitute a significant fraction of
the total turbulence, by the time the dispersion enters the long-time
regime the particles will have moved into a region of the fluid where
the resolvable scale, rather than the SGS, motions dominate the particle
displacements. Proceeding under this assumption, we adopt expressions
for the velocities £ and xn that are consistent with Eqs. (A-16) and (A-9)
and with estimates of the Lagrangian integral time scale T of the SGS
velocity field. The long term statistical behavior of £ and xn that
result from these specifications will be regarded as best estimates of
o2 [cf. Eq. (A-13)].
In the context of the Monte Carlo type of models of diffusion, the
previous comments do not apply because there u' represents the entire
spectrum of turbulent motion. In this case Eq. (A-16) holds for small t
2
and o ~ t applies in the limii
replicates both these regimes.
2
and o ~ t applies in the limit t -> ». The scheme described here
Since the statistical properties of £ and xn are identical, for
simplicity we treat only 5 in the following discussion. A suitable
expression for | is the Markov-1 form
c(tn) = txc
where t = nAt, At is the time step used in the diffusion simulation,
a and 6 are constants, and p (t ) is a random variable with the following
properties:
= 0 (A-18a)
=0 , n f m , (A-18b)
=0 , n f m . (A-18c)
78
-------
The first of these is necessary to satisfy Eq. (A-9).
Through straightforward analyses, we find that with £ given by
Eq. (A-17), its autocorrelation is
and its mean square varies temporally as
?(tn) = c^U^) + B2P2(tn) . (A-20)
Consider for the moment the case of stationary, homogeneous turbulence:
C2(t) = u'2(t) = Constant , (A-21)
—2
where u' is the mean square SGS velocity fluctuation. To ensure that
Eq. (A-21) is satisfied when Eq. (A-17) is used to describe particle
velocities in stationary, homogeneous turbulence, we must set [see
Eq. (A-20)]
a2 + B2 = 1 , (A-22)
p|(t) = u'2 (A-23)
and the initial values of \ must be selected such that
. (A-24)
The Lagrangian integral time scale T of I, which is defined by
79
-------
u' m=0
acquires a particularly simple value in stationary, homogeneous turbulence,
namely [see Eqs. (A-19) and (A-21)],
T, = At am = 7-- (A-26)
^^ I - a
m=0
Thus, in applications to stationary, homogeneous SGS velocity fields with
prescribed values of T and u1 , the parameterization scheme must satisfy
for given At:
a = 1 - f-- , (A-27)
o 1/2
6 = (1 - a ) , (A-28)
and p and e(0) must satisfy Eqs. (A-23) and (A-24).
In many cases of applied interest, the SGS field is stationary but
not homogeneous. In this case,
r(tn) f r(tm) , n t m . (A-29)
If we were to apply the parameterization scheme to problems of this
type, we could continue to use Eqs. (A-24), (A-27), and (A-28) to specify
£(0), a, and g, but we would have to require that p satisfy [see
Eq. (A-20)]
-y
80
-------
2
where u1 (x ) denotes the mean square SGS velocity at the location x
of a particle at time t = mAt on any given realization. Placing this
condition as a constraint on the random variable p is undesirable because,
during the execution of the scheme on the computer, Eq. (A-30) would
require a square root calculation, which is time consuming; and in
regions of highly inhomoyeneous SGS fields, Eq. (A-30) might require
2
(for a given At) that p < 0. For these reasons, we adopt the following
approximation:
= u'2(x)
It can be shown that the error incurred through the use of this approxima-
tion in inhomogeneous turbulence is negligible if
At «
vu'2 II . (A-32)
2 2/3
For example, in the unstable surface layer where w1 ~ (Z/L) [see
Panofsky et al. (1977)], the right side of this expression has a value
of about 100 sec at z = 10 m when l = -25 m and u* = 0.1 m/sec. In short
Eq. (A-17) can be applied to inhomogeneous SGS velocity fields provided
that the time step At used in this equation satisfies Eq. (A-32) where
u'2 denotes the mean square SGS velocity as a function of space.
When the parameterization scheme is implemented, a random number
generator is used to produce values of p . The mean and variance of this
function are specified by Eqs. (A-18) and (A-31), respectively. It
should be kept in mind that the right-hand side of Eq. (A-31) is a known,
local property of the SGS velocity field.
A Gaussian random number generator can be used to produce P£ in the
Monte Carlo type of models applied to neutrally stratified flows because
81
-------
turbulent velocity fluctuations are known to be normally distributed in
this case (Batchelor, 1953). The Gaussian generator should also be appro-
priate in simulating subgrid-scale turbulence in both neutrally and unstably
stratified flows provided that the grid size is much smaller than the
scale of the energy-containing eddies. However, the joint distribu-
tion of turbulence velocities at two sufficiently close points is not
normally distributed under any conditions (Batchelor, 1953). Thus, the
use of jointly normal random number generators in the parameterization
scheme will result in errors in the joint statistics of the particle
pair displacements, but this error is probably no larger than that
incurred from other assumptions in the scheme.
Generating random numbers with a Gaussian distribution requires more
computer time than that needed to generate uniformly distributed random
numbers. Therefore, it is advantageous to use a uniform random number
to index a table of numbers configured in such a way that the numbers
retrieved from the table have a Gaussian distribution with zero mean
and unit variance. This technique requires less than half of the machine
time required to generate Gaussian numbers using conventional algorithms.
In conclusion, l(t ) is given by Eq. (A-17) with a and g
specified by Eqs. (A-27) and (A-28), respectively. The random vari-
able p has a Gaussian distribution with zero mean and variance given
by Eq. (A-31). If the time step At satisfies Eq. (A-32), Eq. (A-17)
can be applied to inhomogeneous SGS velocity fields.
The equivalence of the statistical properties of \. and xn suggests
that an expression similar to Eq. (A-17) should be suitable for producing
xp. In particular, if
82
-------
where a and B are given by Eqs. (A-27) and (A-28), respectively, and
pp(tm) = u'2(xm) , (A-34)
Ao) = So) = u'2(x J , (A-35)
n0
then all of the conditions imposed on xn» namely Eqs. (A-9), (A-10), and
(A-13b), are satisfied. Note that T = T = T, where T is the Lagrangian
X s>
integral time scale of the SGS velocity. The expressions for c and xn
follow immediately from those of £ and x :
n-1
(A'36)
c(tn) =
m=0
m-1
xn(tm) = xn(0) +xn(tk)At . (A-37)
These equations together with Eqs. (A-15), (A-17), and (A-33) constitute
the basis of our parameterization scheme. All that remains to render
these equations operational is to express the functions a and b in terms
of available quantities and to check whether the scheme produces the
proper particle behavior near reflective boundaries.
d. The Velocity Parameterization at Reflective Boundaries
The Markov form of the velocity parameterization we have chosen imparts
a "memory" of approximate length T to the motion of each particle. Conse-
quently, a particle that has a positive velocity, say, just prior to its
striking a reflective boundary, is more likely to have a positive velocity
after impact than a negative one; and hence it is more likely to strike
the boundary a second time than to move away. As a result, the particle
will linger at the boundary. This phenomenon is manifested in diffusion
simulations by fictitiously large particle concentrations at all reflec-
tive surfaces.
83
-------
Actually, turbulent eddies that carry particles into a wall are most
likely dissipated by viscous forces near the wall and, consequently, the
memory they impart to particles is lost. We can approximate this effect
in the parameterization scheme by modifying the expression for £, l»
X , and xn as follows:
if reflection
_ = Pf(tJ > occurred during , (A-38)
5
where p and p are Gaussian random numbers with zero means and variances
._m ^ - - - M ^j~
P^ = p2 = u'2(x ), where u' (x_) denotes the mean square SGS velocity at
the particle's location at time t. Equation (A-38) is used only after
reflection. At all other times, Eqs. (A-17), (A-33), and (A-37) are
used. Tests of this boundary modification have shown that it eliminates
the erroneous accumulation of particles at the boundaries.
e. Specification of the Parameters a and b
When two particles are separated by a distance that is larger than
the largest turbulent eddies, the random motions of those particles are
~~2~ 22 2
uncorrelated. In this case, it is easy to show that I = cr, + 0o (= a
in the present situation). At the other extreme where the particles are
released at virtually the same time and locations,_they never separate
2
(under the action of turbulence alone), and thus, a remains zero for all
time. These two limiting cases indicate that a, as given by Eq. (A-14a),
ranges in value from unity for particles whose motions are totally cor-
related to zero for those whose velocities are statistically independent.
These observations indicate that in the limit of large travel times, as
the mean particle separation becomes large compared with the grid size, the
function a approaches zero. In order to determine the initial value of
a and the manner in which it approaches its final value, we need explicit
P —»j-
expressions for oc and v-.
84
-------
2
Like a, £ is a parameter about which we have little information.
Using Kolmogorov's similarity hypothesis, Batchelor (1950) showed that
if the initial separation £Q of two particles lies in the inertial sub-
range, then for small times:
, (A-39)
where c-, is a constant on the order of 1, G is the energy dissipation rate
of the turbulence, and t is travel tine. In Deardorff's (1974) turbulence
model, A = (AxAyAz) (where Ax, AU, and AZ are the grid cell dimensions)
lies in the inertial subrange; so Eq. (A-39) should apply for all ju less
than A but large enough that viscous effects are unimportant.
Deardorff and Drake (1975) suggest the semi -empirical formula
(MO)
G = 1 + -, - = - , zj.4r , (A-41)
where z is elevation above ground and E is the S6S kinetic energy. By
combining Eqs. (A-14a), (A-16), (A-39) , and (A-40) and by assuming iso-
tropy of the SGS velocity field, i.e., u'2 = 2/3E, we obtain
p/-j/4n\2/3
aQ = a(t0) = 1 - (0.7Gr J(Y) . (A-42)
For convenience, we have chosen c-| = 2/5. It is interesting that the
initial value of a is independent of E. The indicated dependence of aQ
on £Q is consistent with the fact that aQ -»• 1 as ZQ -*- 0, but Eq. (A-42)
cannot be valid for ju » A because it yields large negative values for
a. Although it is possible that a could achieve a small negative value,
in our applications of the parameterization scheme we assume that
0 i a i 1.
85
-------
In addition to deriving Eg. (A-39), Batchelor obtained an expression
for the temporal variation of fc2 during the intermediate time period
2/3 -1/3
t » IQ 6 ' with t small enough that the particle separation is small
compared with the scale of the energy-containing eddies of the turbulence.
We could utilize this information to help describe the behavior of a(t)
between its initial and final values, but this approach would be of very
limited use. Because it is applicable only in the range £Q « £2 « A,
this additional information would be of value only when JL, is several orders
of magnitude smaller than A. Rather than piece together an expression
for a(t), we adopt a single, ad hoc formulation that is consistent with
available information.
Specifically, we note from Eq. (A-42) that the slope of a(t) is zero
2/3 -1/2
at t = 0. Also, we expect that A e must be proportional to the
characteristic time scale of a(t). This information and the known initial
and final values of a(t) can all be taken into account by the following
simple formulas:
a(t) = a exp (-A) , (A-43a)
• /* *'
^0 T2(t')
A = / -?r— , (A-43b)
•(t1) dt1
T(t) = -yr^— , (A-43c)
EV2(t)
where E(t) is the SGS kinetic energy at the centroid of the particle pair
at time t, and X is a constant. Observational studies of mean square con-
centration fluctuations c? in turbulent pipe flows indicate that cVc2
is self-similar over a considerable region downstream of the source.
Lamb (1977) speculates that this is a characteristic feature of dispersion
in decaying turbulence and has shown that it implies that £2/a2 = constant.
[To simulate this condition with our parameterization, we would need
a = constant, which in turn would require that the constant in Eq. (A-43c)
have a very large value. In a future study, we plan to investigate the
value of x in some detail.
86
-------
If the scheme developed here were used in a Monte Carlo model of dif-
fusion in the planetary boundary layer, in place of Eq. (A-40) one would
need the c(z) profile appropriate for boundary layer turbulence (see Wyngaard
and Cote, 1971); and in place of u'2 = 2/3E, one might use the u'2(z)
f
profiles given by Panofsky et al. (1977). This approach is described in
this report in Chapter IV.
f. Recapitulation
For convenience, we summarize below the complete parameterization
scheme for the two-particle subgrid-scale velocities ul and uj>.
n = 1, 2 ;
if a reflection of either particle
occurred during t , <. t < t ;
+ &p?(tm)
otherwise.
if reflection of either particle occurred
during t < t <
"n".'
if a reflection of particle n
occurred during t -j <. t < tm;
otherwise.
if reflection of particle n occurred
during t , <. t < t ;
87
-------
Initial values are c(0) = PC(O), c(0) = 0, *n(0) = Pn(0), and xn(0) = 0.'
Here, Pn(t) (n = 1, 2) and pF(t) are Gaussian random numbers with zero
— — . ^
means and mean square p|(t), and so on equal to the mean square subgrid-
o" ^
scale velocity u1^ at the location of the particle at time t. The param-
eters a and 6 are given by Eqs. (A-27) and (A-28), a and a are given by
Eq. (A-43), and b and 6 are specified by Eq. (A-14b).
3. RESULTS OF TEST CALCULATIONS
In Figure A-l, we show that the scheme gives the correct temporal
variation in the rms separation of particle pairs when we set the param-
eter a(t') = an = constant [see Eq. (A-14a)]. In this case Eq. (A-14a)
~~? 9
gives IL ^ a . This relationship is clearly satisfied by the rms particle
separation values shown in Figure A-l that were determined from the
simulation. Also, by comparing the curves given in Figure A-l for the
case aQ = 0, i.e., uncorrelated motion of the particle pairs, with the
a profiles shown in Figure 6 (Chapter III), one can see that the param-
A
eterization scheme reproduces correctly the limiting relationship
lim &2(t) = 2o2(t)
t -*-"0
The reader is referred to Chapter III for results of single particle
displacement statistics and concentration fields simulated using the
technique developed here.
4. CONCLUSION
We have developed a method of generating pairs of random velocity
fields suitable for simulating particle pair motions in turbulent fluid.
The scheme can be used to simulate the effects of subgrid-scale turbulence
in situations in which one is interested in determining Lagrangian tur-
bulence statistics from instantaneous velocity fields obtained from a
numerical turbulence model; or the scheme can be used as the basis of a
88
-------
CO
100
9
a
7
t
CM
r- __ X
CO O
«*
I
71*1
3 4 S 6 7 8 9 K>
3 4 5 6 7 S 9100
FIGURE A-l. STREAMWISE COMPONENT (a2 - S-Q)l/2 OF THE rms SEPARATION OF
SIMULATED PARTICLE MOTIONS UNDER THE CONDITIONS u = 3.0,
u7^ = 1, AND v7^ - 0.6 (MKS UNITS) FOR TWO VALUES OF THE
PARAMETER aQ
-------
Monte Carlo model of turbulent dispersion, in situations in which one
attempts to determine Lagrangian turbulence statistics from Eulerian
statistics. The inputs to the scheme include the profiles of mean wind
and variances of the turbulence velocities, the Lagrangian integral time
scale, and the energy dissipation rate of the turbulence. With very simple
modifications, the scheme can be used to simulate buoyant particles,
particles with nonzero settling velocities, and particles that are absorbed
or resuspended from boundary surfaces. We are using this scheme in all
of the roles described here.
90
-------
APPENDIX B
USER'S GUIDE TO THE LAGRANGIAN HIGHWAY MODEL
-------
APPENDIX B
USER'S GUIDE TO THE LA6RANGIAN HIGHWAY MODEL
The theoretical aspects of the highway model are discussed in the main
text of this document. This appendix presents a short discussion of the
use of the model. Although the computer program was coded using American
National Standards Institute (ANSI) FORTRAN IV, there are some machine-
dependent portions that the user must check to ensure that the program will
work properly on his computer system. These are discussed in Section 3.
1. DESCRIPTION OF THE INPUT CARDS
Five sets of input cards must be used to run the computer program:
> Overall control cards
> Monitor site cards
> Meteorological data cards
> Source location card
> Miscellaneous data cards.
These sets are discussed in turn below.
a. Overall Control Cards
Four control cards must be placed at the beginning of the input deck.
Card No. 1 contains the number of simulations to be performed. Two param-
eters are entered on this card: the first is the number of the first
experiment (usually 1), and the second is the number of the last experi-
ment. Thus the user wishing to perform three simulations must specify
1 and 3 in Columns 5 and 10 on Card 1. The two variables are entered as
five-digit integer variables. The format for the Card 1 is given in
Table B-l.
92
-------
TABLE B-l. FORMAT FOR CARD 1
Variable Column Format Description
NSTART 1-5 15 The number of the
first simulation
NSTOP 6-10 15 The number of the
last simulation
Card 2, the second card in the input deck, contains the title of the
simulations. The title may be up to 72 characters.
Cards 3 and 4 in the input deck contain 21 variables (12 on the
first and 9 on the second). Each variable and the format for that var-
iable are described in Table B-2.
b. Monitor Site Cards
The calculated PI and P2 kernels (and ultimately the concentration
of pollutant at each monitor site) are stored for each monitor site.
Therefore a simulation is not performed unless there is at least one monitor
site. The maximum number of monitors is 100. The format for the monitor
site cards is given in Table B-3. One card is entered for each monitor
site, and these cards are placed after Card 4. Although the program can
handle up to four lanes, the computer cost increases with the number of lanes.
If the lanes are located symmetrically about the center of the road-
way, the user can create pseudo-stations to reduce the number of lanes to
be treated. (This procedure is discussed in the main text.) The pseudo-
sites will monitor the PI and P2 calculations for the extra lanes.
c. Meteorological Data Cards
Meteorological data are read from two cards. The first, Card 6A,
is described in Table B-4.
93
-------
ZMAX
STRTIM
ENDTIM
TABLE B-2. FORMAT FOR CARDS 3 AND 4
(a) Card 3
Variable Column Format
Description
LFILE
LGRID
LPRNT
NPLT
LMET
NDATA
NX
NZ
XMAX
1
2
3
4
5
6
11-15
16-20
21-25
11
11
.11
11
11
11
15
15
F5.0
Option to save data on an
external file
Option to print calculated
results on a grid
Option to print input data
Option to plot particle
trajectories
Option to print calculated
meteorological data
Flag for comparison of
calculated data and
observed data
Number of cells in the
x-direction
Number of cells in the
z-direction
Length of the region in
26-30
31-35
36-40
F5.0
F5.0
F5.0
the x-direction (in
meters)
Length of the region in
the z-direction (in
meters)
Starting time of the
simulation
Ending time oi
simulation
the
94
-------
TABLE B-2 (Concluded)
(b) Card 4
Variable Column Format
Description
NMON
NSAMP
NITER
1-5
6-10
11-15
15
15
15
NRLMX
TOL
SMPDIM
SEPR
MONSGS
MONWAK
16-20
21-30
31-40
41-50
54
55
15
F10.0
F10.0
F10.0
LI
LI
Number of monitors
Sampling number of conver-
gence test (recommended
value is 7)
Number of 100 iterations
used for sampling (recom-
mended value is 7)
Maximum number of itera-
tions for each simulation
Error tolerance at each
monitor site
Sampling dimensions of
each monitor
Separation between the two
paired particles
Option to activate print-
ing of particle statistics
Option to activate print-
ing of particle statistics
caused by the wake
turbulence
95
-------
TABLE B-3. FORMAT FOR THE MONITOR SITE CARDS (CARDS 5}
Variable
NAMMON(I)
XMON(I)
ZMON(I)
Column Format
Description
1-4
6-15
16-25
A4
F10.2
F10.2
Name of the i-th monitor (if
any
Location of the i-th monitor
in the x-direction
Location of the i-th monitor
in the z-direction
TABLE B-4. FORMAT FOR THE FIRST METEOROLOGICAL DATA CARD (CARD 6A)
Variable
Z0
ZI
ANEHMT
OBUKOV
USTAR
WNDSPD
WNDANG
Col umn
1-5
6-10
11-20
21-30
31-40
41-50
51-60
Format
F5.0
F5.0
F10.1
F10.1
F10.1
F10.1
F10.1
DTDZ
61-70
F10.1
Description
Surface roughness
Inversion height
Anemometer height
Monin-Obukhov length
Friction velocity
Mean wind speed
Mean wind angle perpendicular
to the roadway
Vertical temperature gradient
96
-------
The user may enter experimental data to compare with model predictions
by including Card 6B (see Table B-5) and an external file called TAPE 10.
The meteorological parameters read from the external file will replace
the meteorological parameters read on Card 6A.
d. Source Location Card
The location and description of the roadway is entered on Cards 7,
after the meteorological data. The format for Cards 7 is given in Table B-6.
All distances and heights on Cards 7 are entered in units of meters.
e. Miscellaneous Data Cards
Three cards are required for miscellaneous data, as described in
Table B-7. The data entered on Card 7 control the simulation time.
Card 9 contains parameters that describe the wake turbulence profile.
Card 10 contains descriptions of the vehicles found on each lane.
2. SAMPLE INPUT DECK
A sample input deck is presented in Table B-8. In this example,
the input data are those used in the simulations of the GM experiments.
Users may wish to use these data for their initial simulations.
3. COMPUTER-DEPENDENT PARAMETERS AND SUBROUTINES
The Highway Program requires approximately 44,000 decimal words to
load on the CDC 7600 computer. This number may vary among computer
systems, and some users may find the program too large for their system.
However, certain routines that handle input-output options can be elim-
inated by making them dummy routines. Consider, for example, routine
SGSMON, which is coded as follows:
97
-------
TABLE B-5. FORMAT FOR THE SECOND METEOROLOGICAL DATA CARD (CARD 6B)
Variable
Column Format
Description
(TITLE(I), 1=12,14)
OBUKOV
USTAR
WNDSPD
WNDANG
DTDZ-
4-15
16-25
26-35
36-45
46-55
56-66
3A4
F10.2
F10.2
F10.2
F10.2
Ell. 2
Title of the experiment
Observed Monin-Obukhov
length
Observed friction
velocity
Observed wind speed
Observed mean wind angle
perpendicular to the
roadway
Observed vertical tem-
perature gradient
Next card:
(STREN(I),I=1,4)
Next cards:
(SF6(I),1-1,20)
1-44 4E11.2
Source strength from
each lane of the roadway
1-72 9F8.3 Observed concentrations
(up to nine on each card)
TABLE B-6. FORMAT FOR THE SOURCE LOCATION CARD (CARD 7)
Variable Column Format
Description
NLANE
LAN1LO
DELTAL
TRUCKH
AUTOH
1-10
11-20
110
F10.1
21-30 F10.1
31-40 F10.1
41-50 F10.1
Number of lanes on the
roadway (maximum = 4)
Location of the first lane
(in meters from the
origin)
Distance between the cen-
ters of the lanes
Height of source emis-
sions from a truck
Height of source emissions
from an automobile
98
-------
TABLE B-7. FORMAT FOR THE MISCELLANEOUS DATA CARDS
(CARDS 8, 9, AND 10)
(a) Card 8
Variable
DELTAT
BETA
TIMSCL
A0
Variable
TIMWAK
A0WAKE
BETAWK
UWKCOF
WWKCOF
Variable
VEHWD(I)
VEHSPD(I)
Col limn
1-5
6-10
11-20
21-30
Column
1-5
6-10
11-20
21-30
31-40
Column
1-5
6-10
Format
F5.0
F5.0
F10.1
F10.1
(b)
Format
F5.0
F5.0
F10.0
F10.0
F10.0
(c)
Format
F5.0
F5.0
Description
Time step (in seconds)
Time decay multiplier of A0
Lagrangian integral time scale
Boundary layer initial par-
ticle separation coefficient
Card 9
Description
Lagrangian integral time scale
of the wake component
Boundary layer initial par-
ticle separation coefficient
due to the wake turbulence
Time decay multiplier of A0WAKE
Coefficient of the u-component
of the wake turbulence
Coefficient of the w-component
of the wake turbulence
Card 10*
Description
Width of the vehicle in the i-th
lane
Speed of the vehicle in the i-th
DRGCOF(I) 11-20 F10.1
lane
D!~ag coefficient of the vehicle
in the i-th lane
* One Card 10 is entered for each lane specified in Card 7.
99
-------
TABLE 8-8. SAMPLE INPUT DECK USED IN THE GM SIMULATIONS
1 63
GM EXPERT MEM* WO.
111111 103 30330. 90.1000.1100.
69 7 100 3000 .03 3.0 0.0 FF
11 104.4 9.38
12 104.4 3.31
13 104.4 0.31
21 132.3 9.50
22 132.3 3.63
23 132.3 0.36
31 147.1 9.38
32 147.1 3.03
33 147.1 0.36
41 163.6 9.63
42 163.6 3.61
43 163.6 0.31
51 174.8 9.30
52, 174.8 3.48
33 174.8 0.36
61 189.8 9.60
62 189.8 3.84
63 189.8 0.58
73 209.8 0.36
83 259.8 0.36
IIP 101.0 9.38
12P 101.0 3.31
13P 101.0 0.31
21P 129.1 9.30
22P 129.1 3.63
23P 129.1 0.56
3 IP 143.7 9.38
32P 143.7 3.03
33P 143.7 0.36
41P 160.2 9.63
42P 160.2 3.61
43P 160.2 0.31
31P 171.4 9.30
32P 171.4 3.48
33P 171.4 0.36
6 IP 186.4 9.60
62P 186.4 3.84
63P 186.4 0.38
73P 206.4 0.36
83P 236.4 0.36
41 83.8 9.38
42 85.8 3.51
43 83.8 0.31
444 113.9 9.3
434 113.9 3.6
464 113.9 0.6
32 136.0 9.5
33 136.0 3.3
33 156.0 0.3
60 241.2 0.5
44 82.4 9.58
44 82.4 3.3
43 82.4- 0.3
64 110.3 9.3
63 110.3 3.6
66 110.3 0.3
67 125.1 9.3
68 123.1 3.5
69 123.1 0.3
70 141.6 9.3
100
-------
TABLE B-8 (Concluded)
69
CM EXPERIMENT JO.
111111
69
61
61
63
64
6-3
66
67
68
80
.03
0.8
0.0
2.3
100 30 330. 90
7 100 3000
141.
141.
132.
132.
132.
167.
167.
167.
237.
300.
1
1.0
0.8
22.2
6
6
8
8
8
8
8
3
8
10.3
136, 1
30.
1.0
.3
3.3
0.3
9.3
3.3
0.3
9.3
3.5
0.3
0.3
0
.08
-2.72
Id. 6
0.
.008
3.0 ».0 FT
.20 1.9 2.o2 .0028
1. 1. .8
0.008
101
-------
SUBROUTINE SGSMON
COMMON . . .
RETURN
END
This routine can be made into a dummy by simply eliminating the code
between the second card and the RETURN card:
- SUBROUTINE SGSMON
RETURN
END
If the user is interested only in the predicted results, the follow-
ing routines can be made into dummy routines:
METRRT j
WAKMON / should be kept for testing purposes
SGSMON )
ERPRNT
WRITGD
GRID
CONVT
PLTRAJ
FRAME
AXES
MNSTAT—used only if experimental data are entered
An option of the Highway Program is an off-line plot (CALCOMP) of
the first 20 realizations of the first 100 particle trajectories. If
this option is used, standard CALCOMP routines must be supplied by the
user:
102
-------
NUMBER
LINE
PLOT
SYMBOL
PLOTS
If the user does not have the facility for off-line plotting, he must
create dummy routines for each of the above.
The Highway Program requires the use of a random number generator to
generate a number between 0 and 1. This is done through the use of the
system routine for the random number generator. On the CDC computer the random
number generator function is called RANF(0). This function varies between
computer systems. Therefore, the user must rewrite the statement
RANUM(X)=RANF(X)*GMAX
found in subroutines SGSSET and ADVECT. The card should be changed to
RANUM(X)=(user's random number function)*GMAX
103
-------
REFERENCES
Batchelor, G. K. (1953), The Theory of Homogeneous Turbulence (University
Press, Cambridge, United Kingdom).
(1950), "The Application of the Similarity Theory of Turbulence
to Atmospheric Diffusion," Quart. J. Roy. Meteor. Soc., Vol. 76,
pp. 133-146.
Binkowski, F. S. (1978), "Consistency of a Semi-empirical Theory for Turbu-
lence in the Atmospheric Surface Layer," Atmos. Environ., in press.
Cadle, S. H., et al. (1976), "Results of the GM Sulfate Dispersion
Experiment," GMR-2107, General Motors Corporation, Warren, Michigan..
Chock, D. P. (1977), "General Motors Sulfate Dispersion Experiment:
Assessment of the EPA HIWAY Model," J. Air Poll. Contr. Assoc., Vol. 27,
No. 1, pp. 39-45.
Danard, M. B. (1972), "Numerical Modeling of Carbon Monoxide Concentrations
Near Highways," J. Appl. Meteor., Vol. 11, pp. 947-957.
Deardorff, J. W. (1974), "Three-Dimensional Numerical Study of the Height
and Mean Structure of a Heated Planetary Boundary Layer," Bound. Layer
Meteor., Vol. 7, pp. 81-106.
Deardorff, J. W., and M. A. Drake (1975), "Boundary Layer Data from a Numer-
ical Integration in Three Dimensions," National Center for Atmospheric
Research, Boulder, Colorado.
Deardorff, J. W., and R. L. Peskin (1970), "Lagrangian Statistics from
Numerically Integrated Turbulent Shear Flow," Phys. Fluids. Vol. 13,
pp. 584-595.
Deardorff, J. W., and G. E. Willis (1975), "A Parameterization of Diffusion
into the Mixed Layer," J. Appl. Meteor., Vol. 14, pp. 1451-1458.
Draxler, R. R. (1976), "Determination of Atmospheric Diffusion Parameters,"
Atmos. Environ.. Vol. 10, pp. 99-105.
Egan, B. A., A. Epstein, and M. Keefe (1973), "Development of Procedures To
Simulate Motor Vehicle Pollution Levels," Department of Highways and
Traffic, Washington, D.C. (NTIS PB-233-938).
104
-------
Eskridge, R. E., and K. L. Demerjian (1977), "A Highway Model for the
Advection, Diffusion and Chemical Reaction of Pollutants Released
by Automobiles," Preprint Volume, Joint Conference on Applications
on Air Pollution Meteorology, American Meteorological Society,
29 November-2 December 1977, Salt Lake City, Utah.
Lamb, R. G. (1978a), "Mathematical Principles of Turbulent Diffusion Modeling,"
Lectures presented at the International School of Atmospheric Physics,
13-27 February, 1978, Erice, Sicily, Environmental Protection Agency,
Research Triangle Park, North Carolina, in press.
(1978b), "A Numerical Simulation of Dispersion from an Elevated
Point Source in the Convective Planetary Boundary Layer," Atmos.
Environ., Vol. 12, pp. 1297-1304.
(1977), "Continued Research in Mesoscale Air Pollution Simulation
Modeling--Vol. VI: Further Studies in Modeling Microscale Phenomena,"
EF77-143, Systems Applications, Incorporated, San Rafael, California.
(1976), "The Calculation of Long-Term Atmospheric Pollutant
Concentration Statistics Using the Concept of a Macro-Turbulence,"
Preprint Volume, Third Symposium on Atmospheric Turbulence, Diffusion,
and Air Quality, American Meteorological Society, Raleigh, North
Carolina.
(1973), "Note on the Application of K-Theory to Diffusion
Problems Involving Nonlinear Chemical Reactions," Atmos. Environ.,
Vol. 7, pp. 257-263.
(1971), "Numerical Modeling of Urban Air Pollution," Ph.D. Thesis,
University of California, Los Angeles, California (available from
University Microfilm, Ann Arbor, Michigan).
Lamb, R. G., and J. H. Seinfeld (1973), "Mathematical Modeling of Urban
Air Pollution: General Theory," Environ. Sci. Techno!., Vol. 7,
pp. 253-261.
Lamb, R. G., and VJ. R. Shu (1978), "A Model of Second-Order Chemical
Reactions in Turbulent Fluid," Atmos. Environ., Vol. 12, pp. 1685-1694.
Lamb, R. G., W. H. Chen, and J. H. Seinfeld (1975), "Numerico-Empirical
Analyses of Atmospheric Diffusion Theories," J. Atmos. Sci., Vol. 32,
pp. 1794-1807.
Monin, A. S., and A. M. Yaglom (1971), Statistical Fluid Mechanics
(MIT Press, Cambridge, Massachusetts).
Morse, P. M., and H. Feshback (1953), Methods of Theoretical Physics, p. 795
(McGraw-Hill Book Company, New York, New York).
Panofsky, H. A., et al. (1977), "The Characteristics of Turbulent Velocity
Components in the Surface Layer Under Convective Conditions," Bound.
Layer Meteor., Vol. 11, pp. 355-361.
105
-------
Ragland, K. W. (1973), "Multiple Bay Model for Dispersion of Air Pollutants
from Area Sources," Atmos. Environ., Vol. 7, pp. 1017-1032.
Ragland, K. W., and 0. J. Pierce (1975), "Boundary Layer Model for Air
Pollutant Concentrations Due to Highway Traffic," J. Air Poll. Contr.
Assoc., Vol. 25, No. 1, pp. 48-51.
Riley, J. 0., and G. S. Patterson, Jr. (1974), "Diffusion Experiments with
Numerically Integrated Isotropic Turbulence," Phys. Fluids, Vol. 17,
pp. 292-297.
Shu, W. R., R. G. Lamb, and J. H. Seinfeld (1978), "A Model of Second-Order
Chemical Reactions in Turbulent Fluid—Part II: Application to
Atmospheric Plumes," Atmos. Environ., Vol. 12, pp. 1695-1706.
Taylor, G. I. (1921), "Diffusion by Continuous Movements," Proc. London
Math. Soc., Vol. 20, pp. 196-212.
Teske, M.(E., and W. S. Lewellen (1978), "Prediction of Estimated Highway
Emissions by a Second-Order Closure Model," Draft Report, Environ-
mental Protection Agency, Research Triangle Park, North Carolina.
Wyngaard, J. C., and 0. R. Cote (1971), "The Budgets of Turbulent Kinetic
Energy and Temperature Variance in the Atmospheric Surface Layer,"
J. Atmos. Sci., Vol. 28, pp. 190-201.
Zimmerman, J. R., and R. S. Thompson (1975), "User's Guide for HIWAY, a
Highway Air Pollution Model," EPA-650/4-74-008, Environmental Protection
Agency, Research Triangle Park, North Carolina.
106
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing}
REPORT NO
EPA-60Q/4- 79^023
3. RECIPIENT'S ACCESSION>NO.
ri_5 AND SUBTITLE
A LAGRANGIAN APPROACH TO MODELING AIR POLLUTANT
DISPERSION
Development and Testing in the Vicinity of a Roadway
5. REPORT DATE
April 1979
6. PERFORMING ORGANIZATION CODE
'. AUTHORIS)
8. PERFORMING ORGANIZATION REPORT NO.
R. G. Lamb, H.' Hogo, L. E. Reid
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Systems Applications, Inc.
950 Northgate Drive
San Rafael, CA 94903
10. PROGRAM ELEMENT NO.
1AA815B CA-033 (FY-78)
11. CONTRACT/GRANT NO.
68-02-2733
12. SPONSORING AGENCY NAME AND ADDRESS
Environmental Sciences Research Laboratories - RTP, NC
Office of Research and Development
U. S. Environmental Protection Agency
Research Triangle Park. NC 27711 _.
13. TYPE OF REPORT AND PERIOD COVERED
Final 3/78 - 9/78
14. SPONSORING AGENCY CODE
EPA/600/09
15. SUPPLEMENTARY NOTES
16. ABSTRACT
A microscale roadway dispersion model based on Lagrangian diffusion theory has been
developed and tested. The model incorporates similarity expressions for the mean
wind and turbulence energy in the atmospheric boundary layer, through which the
effects of wind shear and atmospheric stability are taken into account, and a
parameterization of vehicle wake turbulence. Through simple modifications, the model
can be structed to treat particle settling, deposition, and resuspension, as well
as buoyancy of the effluent material. Calm winds, winds parallel to the roadway,
flows around depressed or elevated roadways, shallow mixed layers, and transient of
spatially variable meteorological conditions can all be explicity taken into account
within the framework of the modeling approach.
The model was tested by applying it to the 30-minute experimental periods
reported in the General Motors sulfate study. Of the 1040 predicted values of the
mean concentration of an inert material (SF-), half were found to be within +30
percent of the measured values. The overall correlation coefficient was 0.91. The
computer time (but not core storage) required by the model is directly proportional
to the distance between the farthest receptor and the road. For the studies reported
the model requires on the average 20 seconds of CPU time on the CDC 7600 to
simulate each of the 30-minute General Motors experiments.
17.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.lDENTIFIERS/OPEN ENDED TERMS
c. COSATI Held/Group
*Air pollution
*Roads
*Atmospheric diffusion
*Mathematical models
Lagrangian functions
Development
Tests
13B
04A
12A
14B
13. DISTRIBUTION STATEMENT
19. SECURITY CLASS (This Report/
RELEASE TO PUBLIC
'IED
20. SECURITY CLASS (Thispage)
UNCLASSIFIED
21. NO. OF PAGES
115
22. PRICE
EPA Form 2220-1 (9-73)
107
-------
U.S. ENVIRONMENTAL PROTECTION AGENCY
Office of Research and Development
Environmental Research Information Center
Cincinnati, Ohio 45268
OFFICIAL BUSINESS
PENALTY FOR PRIVATE USE. S3OO
AN EQUAL OPPORTUNITY EMPLOYER
POSTAGE AND FEF.S PAID
US ENVIRONMENTAL PROTECTION AGENCY
EPA-335
If your address is incorrect, please change on the above label
tear off- and return to the above address.
If you do not desire to continue receiving these technical
reports, CHECK Wffffd, tear off label, and return it to the
above address.
EPA-600/4-79-023
------- |