
«»•
^s.
CO
: ^J
^.
"o
III
a
Note: .0 = 2 x 10"8 m"3
2
m = 10 gin/sec
k—5,/mole/sec
2
K,, and K, for source height of 10 m.
k •>. 10
10
10
FIGURE 35. VARIATION OF THE DIMENSIONLESS PARAMETER y AS A FUNCTION OF
RATE CONSTANT k FOR VARIOUS ATMOSPHERIC STABILITIES
CT)

127
where N is the number of "point" sources in the averaging volume. Let A
X
and A denote the separation distances between streets, and assume that
AX » A » w ,
Ay » A » W
y y
Then we have
N = Mk\ + /^\H
/AXAY\/XX +
I WA Vy
Therefore,
A + A
n =
U
WAZA A
x y
The source strength function S that describes the network of streets is
S(r,t) = ^U(r) , (236)
where m is the mass emission rate of each "point" source that constitutes the
network, v is the point source volume given by Eq. (234), and
!1 , if r falls on a street,
(237)
0, otherwise.
The source correlation function G [see Eq. (196)] is therefore given by
U(r' + c) U(r" + c) dc d^ .d^
" ~ ~ x y z
(238)

128
where V = AxAyAz. Making the change of variables
r = r + Ar
and using the fact that G(r',r'+Ar) is independent of r' for all r' lying in
the slab 0< z X X
Therefore,
and Ar
whA_xAy if Ar  ^ nx
v XX
and Ar I = mx
whDV if Ar  = nX
X X
and Ar 
G(r',r'+Ar) = G(Ar) =
A ^— + Dd)(Ar ) d>(Ar )
WX AZ WX AZ vx X; VV y
(240)
if 0 < z'< h , and z' + Az < h ; (241)
6(Ar) = 0 for all z' > h

129
AX
w
Av(r)
~]
FIGURE 36. UNIFORM NETWORK OF STREETS ILLUSTRATING THE CALCULATION
OF THE SOURCE CORRELATION FUNCTION G

130
In this expression,
0 , if Ar  * nx ;
A A
1 , if Ar  = nx and (AT Amx
xx y y
(242)
Similar values apply for 4>(Ar ). Substituting these results into Eq. (195),
we obtain
t) = I ill p(rr',t,t') p(rr'Ar,t,t") G(Ar,t't") dAr dr1 dt1 dt"
p(c+Ar,t,t") G(Ar) dAr d§ dt' dt"
0 0
(243)
P
K>'
+ a'2)
exp
Ar I ,mx
,2
dAry
+ term in (Ar ) + term in cj>(Arx)(Ar• )} dt1 dt" ,
where
**•;?
m2h
(244)
and
= aH(t")
and so forth. Evaluating the expanded term in Eq. (243), we obtain

term in
exp
oi 2 , i2\
2(aH + aH ,
f
2 , ,2x1/2
H + °H } m=
<

132
Determining the mean value Aj is not quite as complex. We find through
straightforward reasoning based on the homogeneity of A, that
Aj(t) =
Total Mass in V
V
mDVt V
AxAya_
Az(t)
jut
wa
X + A
_x y_
A A
x y
(249)
Using this expression and Eq. (248), we can determine the variable i\B [defined
in the present instance by Eq. (189) that is required in the subgridscale
parameterization scheme [Eq. (186)J. Unfortunately, Eq. (248) involves non
elementary integrals that can be evaluated only by numerical integration tech
niques. Such methods will therefore be necessary to implement the present scheme
in an airshed model. For the present, however, it is useful to obtain c.n order
ofmagnitude estimate of the effect of the subgridscale concentration variations
on the chemistry.
According to our scheme, the measure of the importance of the subgridscale
chemistry effect is the degree to which the parameter
1
1AB
A2
AI
differs from unity during the period when the generation of material by sources
and the chemical decay are unbalanced. Through simple, straightforward analysis,
we find that initially
r .1 I
FAB ^ ln~
(250)

133
assuming that A ^ A v, x. Note that rAR is independent of the emission
X jr MD
rate. Thus, assuming an effective emission depth h = 2 m, a street width
w = 10 m, and an average street separation A = 150 m, we find that in a
diffusion model that employs a grid spacing AZ = 20 m in its lowest level,
We know from the work presented in the previous sections that r,, decays
because of turbulent diffusion and ultimately attains a constant value of
unity. The time scale of this decay is given roughly by
X2
where it is assumed that both a and a,, are proportional to t . Since rflR
2_ n MD
is initially much larger than unity, significant subgridscale chemistry ef
fects will arise if T turns out to be comparable to or larger than the time
scale Tp, for example, required for the chemical decay and source emissions
rates to reach an equilibrium. The time scale Tp is equivalent to the time
required for
kA 2 ^ S
where S denotes the source emission rate. In the present problem,
•
S = mD = r
WAAz
and
• 9 9
A2 , 4m V
f\ O O O •
W JTAz
Thus,
/L.M/2
(252)

134
In the present problem, we find that
Tr ^ £v, 400 sec (253)
1 i\i I
o
(using K,, ^ 60m /sec), and
,v _ (254)
E
Subgridscale effects will be most pronounced for fast reactions in which T
7
is small. For the case of the NO0 reaction, we have k ^ 10 £/mole/sec.
Thus, if vehicles emit about 1 gtn of NO per kilometer, on a street network
4
in which 10 vehicles pass each point per day,
TE * 200 sec . (255)
From this estimate and from Eq. (253), we conclude that in areas with street
network densities less than or comparable to that treated here (i.e., streets
150 m apart), and with traffic densities greater than or equal to 10 vehicles
per day, subgridscale concentration variations will have a nonnegligible ef
fect on the simulation of ozone and nitrogen oxide concentrations. To account
for these effects, we must incorporate a scheme such as that developed here
into the simulation model.
3. Isolated Point and Line Sources
The analyses presented in the last two subsections deal with homogeneous
source distributions. The results of those analyses are therefore applicable
only in portions of an urban region where sources are fairly uniformly distrib
uted upwind of a given site. In this section, we briefly consider the magni
tude of subgridscale effects arising from isolated point and line sources, such
as power plants and major highways.
For the purposes of this study, we assume that, in a point source plume,
the mean concentration is described with adequate accuracy by

135
c(x,y,z) = ^— exp  ^exp  ^ , (256)
i . f~r ic • — / I l^.£_E * *
where Q is the emission rate,
2 ,, x
°z  kz  '
a2 = K * ,
y y u
and where the xaxis lies on the plume center! ine with the mean wind parallel
to the xaxis. Assuming that the grid volume V = AxAyAz has dimensions much
larger than the y and z dimensions of the plume, we find after integrating
Eq. (256) that
where d is the effective source diameter. This value of fflR pertains to the
grid cell that contains the source. At a point nAx downstream from the source,
f^g can be found by replacing £n(Ax/d) by £n[nAx/(n + I)AX]. At each point in
space, the value of fAR is constant, because of the constancy of the source
and meteorology, but the value that emitted material "sees" as it moves down
stream decays with time.
In the context of the present problem, the effect of subgridscale concen
tration variations on reaction rate is accounted for in our parameterization
scheme by replacing the true rate constant k by an effective value k' given by
k1 = rABk . (258)
The resulting perturbation in the space averaged concentration due to the subgrid
scale variations is therefore
dc _
~ rr
1 + (kcjt) '

136
where t = AX/U is the residence time of material in each grid cell. With
Eq. (259) reduces to (assuming that kCjt > 1)
AX/TIT
(260)
Typically; AX = Ay, AZ ^ 20 m, and K K ^ 100 n$sec . Recalling that our
analysis assumes that the plume dimensions a and a are much smaller than
AZ and Ay, respectively, we conclude that for sufficiently large u,
* 20 . (261)
This result indicates that the perturbation due to the subgrid scale is very
large for reactions in which (kCjAx/u) > 1. The NO0 reaction in power plant
plumes is one example of a situation that satisfies this criterion. When
kCj(Ax/u) becomes much smaller than unity, the subgridscale effect diminishes.
These results agree with the analyses of randomly distributed point sources
presented in Section IVE1.
Turning finally to the subgridscale effects produced by isolated line
sources, we assumed a line source plume concentration distribution of the form
c(x,y,z) = ^exp (%) (262)
\2oz /
where Q. is the emission rate per unit length of source and
Integrating Eq. (262) and assuming that Az»a2, we obtain

137
1/2,
f  u Az
FAB  l/2
As in the case of the point source, the subgridscaleinduced perturbation
in the concentration will be on the order of
1/2
u AZ
c 1AB A 1/2 K 1/2
z
when
kCj i~ > 1 . (265)
u
Since the line source produces a mean concentration level CT of
Q,
(266)
UAZ
Eq. (265) is satisfied if
2
kQL > . (267)
For the NOCL reaction with k = 10 £/mole/sec, Eq. (267) reduces to
QL > 10"3 gm/sec/m (268)
3
when u = 3 m/sec, Az = 20 m, and Ax = 10 m. If NO emissions are on the order
of 1 gm/km, Eq. (268) is satisfied only by roadways carrying 10 or more vehi
cles per day. Few highways rank in this category. Even for those that do, r.R
is only on the order of unity; consequently, little if any subgridscale chemistry
effect occurs.
We conclude that in contrast to the point sources, line sources of commonly
encountered strengths produce negligible subgridscale chemistry effects, at
least insofar as the N003 reaction is concerned. All .other reactions that are
explicitly treated in air pollution simulations are slower than this one and
accordingly are less affected by subgridscale concentration variations.

V DEVELOPMENT OF A SUBMODEL
FOR RESTORING POINT SPATIAL RESOLUTION
TO GRID MODELS OF URBAN POLLUTION
A. INTRODUCTION
In Chapter I, we have discussed and explained the inability of grid
models (i.e., those in which space and time are discretized) to resolve
features in the concentration field smaller than the discretization interval.
Generally speaking, urban diffusion models of the grid type employ grid
networks with a horizontal mesh size of several kilometers and a vertical
mesh several tens of meters in length. Since most major sources of air
pollution, such as power plants, refineries, and highways, have scales much
smaller than these, there is a great deal of fine structure in the concentra
tion distribution that grid models cannot resolve. Indeed, the locations
and intensities of concentration maxima, the extent and location of the
zone of oxidant depression near roadways, and the pollutant exposures in
street canyons are examples of important pollution phenomena about which
urbanscale grid models can provide no information. Moreover, validation
studies of grid models use monitoring station data that represent time
averages at a fixed point rather than averages over the large volumes re
presented by each grid cell. Consequently, without an estimate of the
amplitude of smallscale concentration variations in the vicinity of each
monitoring station, the validation process can be hindered, because it would
be difficult to ascertain whether discrepancies between the computed and
observed concentrations were due to smallscale spatial variations or to
errors in the model.
These weaknesses have been cited by some in arguments favoring the
development of alternative pollution modeling approaches. Spectral methods
are among the alternatives most often proposed. In this chapter, we develop

139
a "microscale model" (not based on spectral methods) that can be used in con
junction with an urbanscale grid model to achieve point spatial resolution
of the concentration distribution at any specified point. In developing this
microscale model, one fundamental constraint was imposed: The microscale
model must operate totally independently of the urbanscale model. That is,
we sought the capability of examining the nearsource smallscale concentra
tion distribution at any point without having to make multiple runs of the
largescale model. We also wanted the microscale model not to be structured
around or coupled to the urbanscale model. If all of these conditions are
met, the microscale model can be used with any grid model, and the required
calculations can be performed with a minimum of time and effort.
Before proceeding with the mathematical development of the microscale
model, we review, for orientation purposes, some of the basic concepts intro
duced in Chapter I.
Most air pollution models in current use are based on the premise that
the processes governing the concentrations of each pollutant can be described
by an equation of the form
S + R , (269)
. ..
at i ax. ax. ij ax.
I « J
where
u = the ith coordinate component of the mean wind,
K. . = the turbulent diffusivity tensor,
' J
S = the strength of all systematic sources of the pollutant,
R = the time rate of change of the concentration due to
chemical interactions among all of the pollutants present.
In adopting Eq. (269), one tacitly assumes that if all of the parameters in
this equation were known precisely, then the solution c(r,t) of the result
ing equation (and initial and boundary conditions) would correspond exactly
to an errorfree measurement of the concentration at the point (r,t).

140
{Actually, c in Eq. (269) represents a mean value, which we can think of as
a time average over an interval of several minutes, say, depending on the
nature of the diffusivities used [see Lamb (1971)].}
If Eq. (269) is to be solved by numerical methods executed on a digital
computer, this equation must first be "filtered" to remove all smallscale
variations that the grid network cannot resolve. One way to achieve such
filtering is to space average the equation at each point over a volume V
equal to that of the grid mesh. Thus, if the grid network has mesh dimensions
of AX by Ay by AZ, then it can be shown that any space averaged variable de
fined by
x+Ax y+Ay z+Az
(270)
XAX yAy ZAZ
is essentially free of spatial variations unresolvable by the grid network.
[In Eq. (270), r = (x,y,z).] Similar definitions can be written for S and R.
Averaging Eq. (269) in the manner of Eq. (270) and assuming that u and K are
nearly constant on scales comparable to the grid dimensions, we obtain
9C . 77 oC _ a i/ 8C ic in
i u .
J
\
Numerical methods can then be applied to this equation to obtain an approx
imate solution for c.
Thus, a grid model yields the quantity c(r,t) rather than the point
value c(r,t). The difference between these quantities, which we denote by
c"(r,t) = c(r,t)  c(r,t) , (272)
is what we refer to as the subgridscale variations (SSV). Upon subtracting
Eq. (271) from Eq. (269), we obtain the equation governing the SSV:
77 = L' + c u i nil
u k s + R
at i ax. ax.

141
where
S"(r,t) = S(r,t)  S(r;t)
(274)
R"(r»t) = R(r,t)  R(r,t)
Since R is not a function of any spatially variable parameter other than
the concentration, it is clear from Eq. (273) that subgridscale variations
c" can arise only from subgridscale variations S" in the source strength.
Given the usual dimensions of grid meshes used in airshed models (see the
earlier discussion), there are no urban areas in which S", and hence c", is
everywhere identically zero.
It is the mean value c(x,t), rather than the space averaged mean c(x,t)
given by the grid model, that possesses the point spatial resolution and
hence all of the information that we require in air pollution studies.
According to Eq. (272), we can acquire this information from the grid model
by calculating in addition the subgridscale field c". Thus, our approach
to the restoration of point resolution to the grid model is simply to develop
a microscale, or subgridscale, model based on Eq. (273) whose output c"(x,t)
can be combined with that of the grid model c(x,t) to describe the point
field c(x,t) at any arbitrary point.
To illustrate an important property of the subgridscale variation c"
and to give the reader a better physical feel for the nature of c", we
examine this field quantitatively in an elementary example in the next
section.
B. QUANTITATIVE ILLUSTRATION OF THE SUBGRIDSCALE VARIATIONS c"(r,t)
For simplicity, we consider a onedimensional problem governed by the
equation
2
f = K^+ fi(z) 6(t) (275)
9t ^

142
with initial and boundary conditions
c(z,0) = 0
lim c(z,t) = 0
(276a)
(276b)
where K is the turbulent diffusivity, assumed to be a constant, and 6 is the
delta function. The solution of Eqs. (275)and (276) is
c(z,t) = (4TrKt)~1/2 exp fz
4Kt
(277)
Now let
z+Az
c(z,t) =
2Az
c(z',t) dz1
ZAZ
Averaging Eq. (277) in this manner, we obtain
c(z,t) =
4AZ
erf
 erf
where
(278)
(279)
erf(x) = 2= f
•v7r J
X
dx
0
is the standard error function. Note that Eq. (279) is the solution of
the equation obtained by averaging Eq. (275); i.e.,
(280)

143
where
I 1 , AZ £ z < AZ
u(z) =
I 0 , otherwise
Let us assume that Eq. (280) is a grid model representation of Eq. (275)
and that c, as given by Eq. (279), is the model's output. The subgrid
scale variations are therefore [from Eqs.. (277) and (279)]
c"(z,t) = = exp  f L „ erf
4AZ 4K
(281)
We can see from this expression that the amplitude of the SSV decreases as
the discretization interval AZ in the grid model is made smaller. More
over, for fixed AZ, the SSV in this example decreases with time. The latter
effect is shown graphically in Figure 37. The figure also shows c for com
p
parison. It can be seen that for the earlier travel time (i.e., t = 0.031AZ /K) ,
c" is up to three times as large as c and is therefore not a negligible
variation.
To obtain a more concise description of the relative magnitude of the
subgridscale variation c", we note first that c" is largest at the source
location (i.e., z = 0). Thus, the maximum amplitude of c" relative to c
can be expressed by
• •
From the expressions for c" and c, we obtain
p(t)  j=. [a* erf (^)j "'
where
(283)
\/4Kt
* = "

144
1.8
1.7
1.6
1.5
1.4
1.3
1.2
1.1
1.0
0.9
0.8
c
o
S 0.7
J
+J
§ 0.6
c
o
LJ
0.5
0.4
0.3
0.2
0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
A
="f
C (a* = 0.25}
\
\ 0.4 \ oTe oTs uo
\ N.
\ ^\
\ ^
\ /
\
1.2
Amplitude
FIGURE 37. ILLUSTRATION OF THE SUBGRIDSCALE CONCENTRATION VARIATIONS
ARISING FROM AN INSTANTANEOUS POINT SOURCE

145
That is, o* is the approximate halfwidth of the pollutant cloud at time t
normalized by the grid mesh size Az. It is evident from Eq. (283) that a*
is the key parameter that determines the relative size of the SSV in this
instance. Upon evaluating Eq. (283) for several values of a*, we obtain
p = .10 when a* = 0.1, p = 0.34 when a* = 1 , and p = 0.08 when a* = 2. From
these values, we conclude that, in the case of a point source, the SSV can
be neglected relative to c after the pollutant cloud has grown to a width
of about four grid intervals.
These results take on added meaning when one realizes that the one
dimensional equation considered here is an approximate model of a two
dimensional steadystate plume. In this case, time is related to the down
wind distance x from the plume source by the relationship
x = ITt , (285)
where u is the mean wind speed (directed along the xaxis). Translated
into terms of the twodimensional plume, the conclusions reached above re
garding the relative magnitude of the SSV indicate that SSV effects will
be significant within a distance of approximately
x = j (286)
from the plume source. In the grid model, this distance x will represent
e* (287)
grid intervals from the source. Thus, for fixed u" and K, the number of
grid cells affected by SSV actually increases as the grid size Az is in
creased, albeit both c" and c become insignificant sufficiently far from
the source.
Although we considered only the maximum value of c" (i.e., c"(0,t) .in
the analysis above, we can expect (see Figure 37) that the SSV will become
negligible everywhere in space on the same time scale as that of the
maximum.

146
C. DERIVATION OF A SUBMODEL OF c" FOR LINEARLY REACTIVE POLLUTANTS
1. Discussion of General Operational Problems
As stated in the introduction to this chapter, we seek to obtain c"
from its governing equation [Eq. (273)] and to use this result along with
that given by the grid model to obtain the point mean field c(x,t). Since
our earlier discussion of the equation governing c", i.e., Eq. (273), was
very general, it did not reveal several of the operational problems that
arise in solving this equation in given situations. We discuss these below,
within the context of a simple problem, before proceeding with the develop
ment of specific models.
To keep the mathematics as simple as possible, we consider here only
a single pollutant that decays at a nonlinear rate. The equation governing
the concentration of this pollutant in the urban atmosphere is assumed to be
H+ LC = kc2  k + S , (288)
where [_ is the operator
L s «1 SJq ^ KIJ alj • <289>
u. is the mean wind speed, K• is the turbulent diffusivity tensor, and
/ 2
(c1 ) is the turbulent concentration fluctuation term discussed in Chapter III,
Although this example is mathematically much simpler than that for photo
chemical air pollution, which involves a system of coupled nonlinear equa
tions, the problem posed by Eq. (288) has many of the essential character
istics of the more complicated real situation. Thus, techniques developed
for a pollutant that undergoes a secondorder reaction can be applied to
more general situations through straightforward extensions.
Suppose that a grid model for the pollutant cited above is to be
developed. As outlined earlier, the procedure for deriving the working
equation of the model is to perform at each point a spatial average of

147
Eq. (288) over a volume V equal to that of the grid cells. In this way,
we obtain the working equation
*"» '—»
f + [_c = kc2  kc"2  k + S , (290)
d U
where c and S are spatial averages and c" is the subgridscale concentration
variation. Note that it is the spatial average mean square turbulent fluc
tuation term <(c )> that enters into the grid model equation. Upon subtract
ing Eq. (290) from Eq. (288), we obtain the equation governing c":
9t
:"2  2kcc" + kc"2  k" + S" , (291)
2
where " represents the subgridscale portion of the mean square fluc
o
tuation field (c1 >. The second term on the right side of Eq. (291) represents
the chemical interaction of the subgridand supergridscale concentration
fields. Since the latter is a known quantity at each point (it is given by
the grid modol), Eq. (291) contains only one dependent variable, c" (after
t^^>
c"^ and (c'^) have been parameterized), and this equation can thus in
principle be solved.
The nonlinearity of Eq. (291) dictates that its solution can be approxi
mated only by using numerical integration techniques. It is to be expected,
however, that the implementation of these techniques will require the solu
tion of operational problems unlike those associated with the modeling of c.
This requirement arises because point and line type sources give rise in S"
to delta functions, which are not amenable to treatment using discrete grid
numerical methods. We could, of course, attempt to evaluate Eq. (291) on
a fine grid, but this would only create new subgridscale problems. To
circumvent these problems, we develop a new technique later in this chapter.
In contrast, in cases where the pollutant is linearly reactive, the
implementation of the microscale model is quite simple because the governing
equation, Eq. (291), reduces to the much simpler form
3r "
+ Lc" = kc" + S" , (292)

148
and if the pollutant is inert, the equation reduces further to
^+ Lc" = S" . (293)
In the remainder of this section, v/e develop the known analytic solu
tion of Eq. (292) into a working model of c" that can be used with any grid
model of CO, S0?, or other essentially linear pollutants.
2. Derivation of the Working Equations for a Microscale Model of
Linearly Reactive Pollutants
Consider a pollutant that decays linearly at a rate described by
^ = kc , (294)
where k is the rate constant. (Note that inert pollutants are represented
by k = 0.) The equation governing c" in this case is Eq. (292), the solu
tion of which is of the form [see Lamb and Neiburger (1971)]
t
c"(r,t) = / / p(r,tr't') exp ~k(t  t1)] S"(r',t') dt1 dr1
'O
+ f p(r,tr't0) exp [k(t  tQ)] c"(r',t0) dr'
(295)
where p denotes the Green's function of Eq. (292) and its boundary conditions,
For simplicity, we assume that c"(r,tQ) = 0. in which case the last term in
Eq. (295) vanishes. Our aim now is to derive formulas for c" pertinent to
general source configurations S" in conditions of open terrain where the
form of p is rather simple. Consider first the source function S".

149
The source distribution S can be expressed in the general form
Spa(t) 6(^  rJ + 2 SL3(t) 6(Z  Z3} 6[y  S3(X)
P (296)
where S denotes the emission rate (mass/time) of the ath point source (of
pa
which N are present), r is the position of the ath source, and 6 is the
delta function. The line sources represented by the last term in Eq. (296)
have strengths S\_ per unit length and lie along the curves y  Sg(x) at
elevation Zp. From Eqs. (296), (274), and (270), we obtain
N
S"(r,t) = S (t) 6(r  ra) + SLg(t) 6(z ' V
Q. ..
T ~} V
t J IS. it / \ / OQ "7 \
AxAyAz ijk ~ ' ^ '
where Q. ., is the total mass of pollutant emitted per unit time in the grid
1 J K
cell (i,j,k) and
!1, if r lies in the grid cell (i,j,k);
(298)
0, otherwise.
Upon substituting Eq. (297) into Eq. (295), we obtain an expression for c" of
the form
N M L
c"(r,t) = r. P (r,t) + > L (r,t)  > V (r,t) , (299)
' a ~ .«•! * p ~ i •• Y ~
a=l 3=1 Y=1
where P is the contribution to c" from the ath point source, L. is the con
a p
tribution to c" from the 3th line source, and where V is the collective con
Y
tribution of all sources in the yth grid cell. Below we derive the functional
form of each of these terms, assuming that in the open terrain and quasisteady
flow regimes of interest to us here the kernel p in Eq. (295) has the Gaussian
form

150
p(r,tr',t') =
where
exp 
(x
 X1  UT)
9 2
2o
X
2 (y
 y' 
9 2
2ay
v,)2]
exp
(2  Z')2
9 2
z
+ exp
(z + z1)2]
9 2
2a
T = t  f,
° = F(T)»
(300)
(301a)
(301b)
u, v = constants,
w = 0.
(301c)
(301d)
This form of p contains the implicit assumptions that the earth acts as a
reflective barrier to the pollutant particles, that the turbulence is homo
geneous, and that no lowlevel inversions exist. Assumptions (SOlc) and
(SOld) are not generally valid, but are acceptable approximations over the
relative small areas, comparable to a grid cell, in which c" is finite. To
maximize the accuracy of Assumption (301b), we plan to use the "optimal"
a (t) and a (t) profiles developed in Chapter II.
X Z
a. The Point Source Contribution, P
Using the delta function's property
/ f(x) 6(x  XQ) dx = f(xQ)
a

151
/•*
«(r'V = Spa(t)J Pxy(x'y'WT)
'T) dT '
(302)
where
Pxy(x,y,Xa,ya,T) = [27rax(T) ay(i)
exp
o
(x  Xa  UT)
2a2(T)
v '
(y  ya 
2a2(i)
v
2
(303)
and
pz(z,za,i) = /2T
1
exp
r ( }
az
+ exp
(Z + Za) ]
2a2(r) J
(304)
We should mention that the assumptions of a reflective earth and the absence
of an upper level inversion can be relaxed by substituting the appropriate
expression for p in Eq. (302). Equation (302) is the first of the set that
we need to model c".
b. Line Source Terms
We first assume that all line source segments are straight lines lying at
ground level (i.e., z = 0, 3 = 1,2, ..., M) and that the strength of each seg
P
ment is constant along its length. Suppose that the 3th line source has
length I that it is centered at (xn, yj, and that it makes an angle 0_ with
p P P p
the yaxis as shown in Figure 38.
FIGURE 38. PARAMETERS REQUIRED TO SPECIFY A FINITE LINE SOURCE

152
The line source term L which represents this source in the c" formula
[Eq. (299)J, can be shown to be [see Lamb and Neiburger (1971)]
0. _ C
r ( [(x  UT  x ) cos 0R+ (y  VT  y ) sin 0 ]
Mr.*) = Slo(t) / exp {
2a
1
erf
(y  VT  y ) co.s 0  (x  UT  x«) sin e +
P P P p
 erf
• K
(y  VT  y ) cos 0  (x  UT  xj sin 0Q  *•
PZ(Z,0ST)
(305)
where a = ax = ay and where pz is given by Eq. (304), with
the general line source term needed to compute c".
= 0. This is
c. Volume Source Terms
The volume source terms V arise from S, each volume source having the
dimensions Ax • Ay • AZ of the grid cells in the urbanscale model. It is
not difficult to show that the volume source term V is given by
V (r t) = ^
Y^' ^ 8AV
^ X  UT + ^\ /X  X  UT
1 ^j _ erf ' ^
 VT
erf I
/2
 erf
fy  y  VT 
/2
Ferf(A)  erf(B) + erf(D)  erf(c)l
(306)

153
where
z  z
/2 o /2 a
(307)
Zj _____ 747 4
1  O i. ~ i. ~ ~7T
B =
and (x ,y , z ) denotes the center of the yth grid cell. Equation (306) is
the volume source term required in Eq. (299) to calculate c".
3. Design of the c" Model Computer Program!nteraction
with the Source Inventory Data Bank
We have derived the basic equations required to calculate c". The next
step is to put these results into the form of a subroutine that the airshed
model can call to compute concentrations c at any given point. As previously
stated, the c" model must not interfere with the operation of the full urban
scale model, and it must not require the latter to perform special calculations
In short, the c" subprogram must be capable of calculating c" given only the
following information (which is always available in the grid model itself):
(u, v) = horizontal components of the mean wind in the
vicinity of the receptor,
z./L = local stability parameter,
AX, Ay, AZ.= dimensions of the grid cells,
u = local friction velocity,
z. = mixing depth in the vicinity of the receptor,
*c = value computed by the grid model for the
receptor point,
(x,y,z) = coordinates of the receptor.
Thus, if MICRO is the name of the c" model, the call for it by the grid model
might look something like the following:

154
CALL MICRO (X,Y,Z,T,ZlOVL,USTAR,ZI,UBAR,VBAR,r
DELTAX,DELTAY,nELTAX,CTILDE,CPP)
where CPP = c"(x,y,z,T).
The MICRO program must also have access to the source inventory, where
each point and line source segment, or link, is stored separately. In partic
ular, the following source data are needed to compute c":
XPS(I)j
YPS(I)> = (x,y,z) coordinates of Ith point source.
ZPS(I))
XLS(J) ) _ (x,y) coordinates of the center of the Jth line
YLS(J)  line source (ZLS = 0 assumed).
THETA(J) = angle of the Jth line source (see Figure 38).
XLNGTH(J) = length of the Jth line source.
SP(I) = strength (mass/time) of the Ith point source.
SL(J) = strength (mass/length/time) of the Jth line source.
Using all of the above information, MICRO computes c" in a manner described
roughly by the flow chart shown in Figure 39. We have completed the steps
labeled SIGMA and CALC in the flow chart; Appendix C presents listings of
their source programs. Steps SEARCH, SORT, and DOMAIN have not been completed,
and without these modules, the MICRO subprogram lacks the capability of search
ing the source inventory itself. However, if the user provides the relevant
source information, the existing steps SIGMA and CALC can compute c" at any
given point.
In its present form, SIGMA consists of two FORTRAN functions, SIGMAX(T)
and SIGMAZ(T), which are used by the subprogram CALC. The SIGMA functions
calculate a and a for a given travel time T and for given stability, mixing
X £.
depth, and friction velocity conditions (which are passed through common blocks
to the function programs). To evaluate the integrals entering in the equations
derived earlier for c", the CALC subprogram uses a Monte Carlo quadrature
technique, described theoretically in Appendix D. This technique allows infi
nite flexibility in the ability of the CALC program to evaluate c" both near

155
CENTER J
Using given values of USTAR, ZI0VL, and II,
convert the stored nondimensional universal
functions OX(T ), ay(r*), and a^t*) into
the dimensional forms ox, ay, and az.
\/
Using ox, ay, az, UBAR, VBAR, X5 Y, and Z,
compute the boundaries of the "domain of
influence," i.e., the region containing all
sources that might affect c"(X,Y,Z,T).
Search the source inventory for all point
and line sources that lie in the domain of
influence.
\/
Store location and strength data for all
sources that may have a significant effect
on c"(X,Y,Z,T); reject all other sources.
\/
Calculate c"(X,Y,Z,T) using Eqs. (299),
(302), (305), (306), and the source data
listed by the SORT step.
f RETURN J
SIGMA
DOMAIN
SEARCH
SORT
CALC
FIGURE 39. FLOW DIAGRAM' OF THE c"MODEL ROUTINE MICRO

156
and far from point and line sources. It also requires minimal complexity
in programming and permits control over the accuracy of the integral eval
uation.
Test calculations performed with the MICRO model in its present form
for the problem portrayed in Figure (40a) are presented in Figures (40b)
and (40c). Although these figures are selfexplanatory, we wish to call
attention to the size of c" in the immediate vicinity of the line source
in this calculation: It is over two orders of magnitude larger than c at
the same point [cf. Figures (40b) and (40c)J.
One of the main problems that we anticipate in the implementation of
microscale models such as MICRO into fullscale urban diffusion models is
the incompatability of the source inventory data bases that the two kinds
of models require. That is, in grid models, such as that developed by SAI,
the source inventory consists of total emissions from each grid cell rather
than emissions from each point and line source within a cell. The MICRO
model requires data in the latter form, which is not derivable from the
existing source inventory.
A logical solution to this problem is to collect raw data in the form
of point and line source emissions and to store it in this form. It is then
a simple matter to compute total emissions from any grid cell network for
use in a grid model. This approach was employed in a pollution modeling
study of the Los Angeles basin conducted by Lamb and Neiburger (unpublished),
Data from this study consist of the coordinates and strengths of each of ap
proximately 6000 roadway links in the Los Angeles area for the year 1966.
These data are stored on discs and can be converted, using simple programs,
into emissions from any desired grid cell network. We have procured this
data bank and will use it in the final developmental phases of the MICRO
program. Our chief interest is in using this data to test the DOMAIN,
SEARCH, and SORT modules, which will allow the MICRO program to interact
with the data bank for the purposes of computing c" at any desired point.

u = (2,3)
* SLIGHTLY
UNSTABLE
CONDITIONS
MEAN WIND VECTOR
10METER
STACK
.GROUNDLEVEL
LINE SOURCE
•GRID CELL
(200 x 200 x 20)
(a) Source Configuration Showing the
Conditions of a Test Problem for
Solution Using the CALC Routine
Listed in Appendix C.
0.02
0.07
0.01
(b)
c" (x 10J) Computed by CALC for
the Problem Given in Part (a).
Note the large magnitude of c"
near the line source. The
elevation of the point source
results in positive c" only at
large distances.
(c)
c (x 10 ) as a Perfect
Grid Model Would Compute
It for the Problem Shown
in Part (a).
FIGURE 40. THE CONTOURS OF c" ARISING FROM A SINGLE POINT AND LINE SOURCE, AS COMPUTED BY CALC
en
—i

158
The final version of MICRO will enable the user of a grid model to restore
pointscale spatial resolution to his diffusion model at any desired point.
The following sections describe the DOMAIN, SEARCH, and SORT modules that
we are presently developing.
a. DOMAIN
We demonstrated analytically in Section VB that c" decays'with travel
time, even in the absence of chemical decay. Let the characteristic time
scale of this decay be T. Thus, sources located a distance
I > [02 + v2]1/2 T (308)
upwind of the receptor point (XR,YR,ZR) have a negligible effect on c". Fur
thermore, because of the limited horizontal and vertical dispersion rates, all
sources downwind of (XR,YR,ZR) and those outside a wedgeshaped sector upwind
of (XR,YR,ZR) have negligible effect on c". By ignoring all such sources, the
MICRO program can optimize its calculation speed.
We showed in Section VB that the decay time scale T of c" is approximately
equal to the time required for a point source cloud to grow to the size of a
grid cell. Thus, if <*x = axtbx> °. = avtby» and az = aztbz> we have
(309)
)
The upwind dimension £ of the domain of influence is therefore
L = VT , (310)
and T is given by Eq. (310).
The width of the domain of influence is proportional to a (t), assuming
that a  a . We conclude that the domain of influence is the wedgeshaped
x y
region shown in Figure 41.

159
(XR,YR :
DOMAIN OF
INFLUENCE
GRID CELLS
> x
(a) Horizontal View of the Domain of Influence
•AZ
RECEPTOR
= UT
(b) Vertical Cross Section View AA of the Domain of Influence
FIGURE 41. HORIZO.NTAL AND VERTICAL CROSS SECTIONS OF THE
DOMAIN OF INFLUENCE ON c"(XR,YR)

160
b. SEARCH
Once the boundaries of the domain of influence have been established,
the next step is to find and list all point, line, and volume sources that
lie within it. As indicated earlier, the coordinates of the ith point
source are stored in the data bank as XPS(I), UPS(I), ZPS(I). Four param
eters are stored for each line source. By testing the coordinates of each
source, we can determine whether that source lies wholly or partially within
the domain of influence.
The speed of the point and linesource searching operation can be opti
mized if sources in the data bank are catalogued according to the grid cell
in which they lie. If this is the case, only those blocks of data filed in
locations corresponding to grid cells within the zone of influence need be
examined. Determining the volume sources within the domain of influence 12
a simple matter, since these correspond to the grid cells. Note that even
if only the edge of a volume source (grid cell) falls within the domain of
influence, the entire source must be treated. The extent to which that source
affects c" at the receptor is automatically accounted for fin the c" equation.
c. SORT
The efficiency of the c" calculation can be increased by considering only
those sources that together have the largest influence on c" at the receptor.
For example, the SEARCH step may find a total of 10 point and line sources
within the domain of influence, but only three of these together may be re
sponsible for 90 percent of the observed magnitude of c". In this case, neg
lecting the other seven sources would cut computation time by 70 percent at
the cost of only a 10 percent error in c".
Because the SORT process itself consumes computing time, the usefulness of
this step is limited to those cases where there are many point and line sources
within the domain of influence. To determine which sources to neglect, we need
an algorithm that permits rapid estimates of the approximate effect on c" of a
given source. The following is such an algorithm.

161
To each point and line source, we assign a number RANK. For the ith
point source, RANK is given by
(311)
where
d  [(XPS(I)  XR)2 + (YPS(I)  YR)211/2 . ('31 2)
For the j^th line source,
RflNK = _.  ( ,
" +
where
THETA(J) + TAN' . (315)
d = [(XLS(J)  XP)2 + (YLS(J)  YP)2] , (314)
w = XLNGTH(J)*SIN
With a value of RANK assigned to each point and line source, we now form'
the sum
RANKT  SRANK , (316)
where the summation is over all point and line sources and where we arrange
the sources in order of decreasing values of RANK, with point and line sources
mixed as their values of RANK dictate.
In the CALC step, the sources are treated in order of their RANK values,
the largest being treated first. As we proceed through the list of sources,
we keep a running sum of the RANK values of the sources treated. Let this sum
be denoted by

162
K
RUNSUM(K) = ]T RANK(k) , (317)
k  1
where K represents the number of sources out of the total of KMAX point and
line sources in the domain of influence that have been treated. Note that,
by definition, RUNSUM (KMAX) = RANKT. Thus, when
RUNSUM(K) ,.1R>
RANKT  B (6lti)
where 0<3<1, we terminate the calculation of the point and line sources con
tributions to c" and neglect all sources with labels between K and KMAX. The
value assigned to 3 depends largely on the estimated overall accuracy of the
model; it is generally unnecessary to set 3 = 1 because errors in the model
itself do not warrant the added computing time required.
A printout of the list of the sources actually treated and their RANK
values will also be of aid in assessing the polluting potential of each
source in the vicinity of the receptor.
D. DEVELOPMENT OF MATHEMATICAL METHODS FOR MODELING SUBGRIDSCALE
CONCENTRATION VARIATIONS c" OF NONLINEAR POLLUTANTS
In the last section, we derived solutions of Eq. (292), which governs
the c" distribution of linear pollutants. Needless to say, obtaining solu
tions of the more general equation [Eq. (291)], which governs nonlinear species,
is a much more difficult process. And modeling the complete system of photo
chemical pollutants is even more complex. We have only begun to consider the
nonlinear problem. We outline here a method that we have developed and are
considering for constructing a microscale model of nonlinear pollutant species.
We demonstrate the method in the context of Eq. (291) because it is mathemati
cally much simpler than the equations describing the complete, kinetic mechanism,
but yet it retains the essential nonlinear property of those equations.

163
The obvious approach to solving Eq. (291) is to express the equation
in finite difference form and to integrate the resulting expression numer
ically. However, this method only creates a new subgridscale problem:
because due to the presence in S" of delta functions arising from point and
line sources, subgridscale concentration variations will exist in any dis
crete analogue of Eq. (291). Even if it were possible to achieve a grid
network fine enough to resolve the major variations in the concentration
field, operational problems would still remain. For example, because of
the arbitrariness of the locations and orientations of point and line
sources, a variable grid network would be required to optimize the accuracy
of the numerical integration process. But this might lead to difficulties
with computational stability, boundary conditions, truncation error, and the
like, and would certainly complicate the computer program.
To circumvent all of the problems, we develop in this section a rela
tively new type of pollution modeling equation. This model is not based on
any radically new theory. Rather, it is based on an equation that is an
intermediate step in a sequence of mathematical operations that leads from
the general Lagrangian equation of turbulent diffusion Eq. (16), through the
wellknown Gaussian puff and plume formulas, to the classical diffusion equa
tion [Eq. (291)J. One of the unique features of the new equation that enables
it to avoid the problems listed above is that it is continuous in space (like
the diffusion equation), but discrete in time. In the remainder of this sec
tion, we derive the new model, briefly discuss its place in the scheme of
existing modeling equations, and outline a set of ways in which this model
can be implemented in the calculation of c" of nonlinear pollutants. The
actual implementation process will be undertaken in a later study.
1. Derivation of a DiscreteTime, ContinuousSpace Diffusion Equation
The starting point of the derivation is the Lagrangian form of the dif
fusion equation introduced in Chapter II, i.e. ,

164
= ff p(r,tr',t') S'(r',t') dt1 dr' + /p(r,t r' ,tQ) dr'
0 (319)
where S is the source strength function. In cases where linear chemical decay
occurs at a rate
a£=k(t)c , (320)
where k(t) is the rate "constant," it can be shown [see, for example, Lamb (1971)]
that the proper form of the kernel p entering in Eq. (319) is
r t i
p(r,tr'Jt1) = p'(r,t r',f) exp  I k(t") dt"
LJt, J
(321)
where p1 is the probability density that the fluid will advect a particle
from (r',t') to (r,t). Note that p1 is a property of the fluid motions alone
and is independent of the chemistry. Through a series of assumptions and math
ematical operations, it can be shown that the Gaussian plume formula, the puff
model, and the classical diffusion equation are all derivable from Eq. (319).
It is in the derivation of the diffusion equation that the discretetime,
continuousspace model arises.
Thus, let the time axis be divided into equal intervals of length At, and
t denote the time t  nAt. For br
the derivation by writing, for example,
let t denote the time t  nAt. For brevity, we omit the reference to time in
Consider now the conditional probability density
p(rnlrn_r rn_2,...,r0) . "(322)

165
that a particle is at r at time t given that it was at r , at time t , ,
~n n  ~n 1 n 1
at r ? at time t _?, and so forth. If an interval At exists such that
ML. lie.
p(r r ,,r 0,...,rn) = p(r r ,) , (323)
p ~n'~nl ~n2' '~0 K ~n ~nl ' v '
then the particle motions are said to consitute a Markov process in the dis
crete time frame because any stochastic process whose probability density
satisfies a relationship of the form of Eq. (323) is called a Markov process.
An important property of a Markov sequence is that its probability density
p(r r ) satisfies the ChapmanKolmogoroff equation
~n ~m
p(rjrnl) =/p(rnlr£) p(r£ rj dr£ . (324)
where n>£>m are any integers.
The physical significance of the Markov process in the context of a tur
bulent diffusion problem warrants comment. Suppose that the turbulent eddies
that are responsible for the random motions of the pollutant particles have a
time scale T. A reasonable estimate for T is the Lagrangian integral time
scale of the turbulence defined by
(t) v (t + s) d£ ' (325)
0
where v'(t) is the velocity component at time t of a fluid particle. (If the
turbulence is stationary, then T. is independent of the time t.) The time
scale T. can be regarded as a rough measure of the lifetime of the socalled
energycontaining eddies.
Suppose now that a particle is released in the turbulence and that its
velocity is recorded at intervals At = T. /10 apart. We can expect that if
the velocity is, say, positive at times t , t +, , and t +?, then there is a
greater chance of observing a positive velocity at time t + ~ than a negative
velocity; and conversely, if the velocity is negative at the first three times,

166
then it is more likely to be negative than positive at time to In con
trast, if the particle velocity is recorded at intervals At = 10 . apart,
then we can expect the velocity at time t + ~ to be independent of that at
t ? and all previous times. The idea here is simply that a particle has
a "memory" of some length T in the sense that the behavior of the particle
at any time t is seemingly influenced by its history during the interval T
just prior to the time t. This apparent memory is imparted to the particles
by the organized or systematic behavior that each eddy exhibits during its
lifetime. It is the randomness of the birth and death processes of the
eddies and of the transferral of a particle from one eddy to another that
causes a particle ultimately to "forget" its distant past history. This
type of reasoning leads us to expect that the probability density of the
particle positions will satisfy the Markov condition [Eq. (323)] if the
discretization time interval At is chosen such that
At » T. , (326)
where T, is given by Eq. (325).
,
For reference, we rewrite here the general equation [Eq. (319)] using
the notation adopted above:
r , r/"'
=/p(rrl r0) dr0 +11 p(r,lr'> S(r',t')dt' dr'
^
p(r r') S(r',t') dt' dr
(327)
The reason for the split integral in this equation will become clear later.
For any time t1 < t ,, we have from the Markov assumption
n~ i
p(rnlrn_,r') =p(rJv,) •

167
Thus, from the general relationship
f(X X ) =
n m
X£'Xm>
(329)
where X , X,,, and X are any three random variables with probability density
11 XL 111
functions f, it follows from Eq. (328) that
p(rn r1) =
(330)
Thus, p(r r') satisfies a ChapmanKolmogoroff type of equation, even through
r' can be the particle position anywhere on the continuous time axis prior to
Vr
Using Eq. (324), we can express the first term in Eq. (327) in the form
r0>
(331)
Similarly, Eq. (330) permits the second term in Eq. (327) to be written as
jj" p(rnlr') s(r',t') dt dr = f p(rnlrn_,)
L/N
//nl
' P(rn_, r') S(r',t') dt' dr1
^ tn
(332)
Combining Eqs. (331) and (332) and comparing the result with Eq. (327) for
/c(r,t ,)") , we conclude finally that, by virtue of the Markov property,
\ ~ n1 '
Eq. (327) is equivalent to

168
^p(r,tnr',tn_1)dr'
P(r,tn r',t') S(r',t') dt1 dr' . (333)
'n1
This is the discretetime, continuousspace diffusion equation.
From the modeling standpoint, Eq. (333) has several very useful assets.
The first is that it allows concentrations at time t to be calculated from
those at time tAt and from the source emissions that have occurred during
the interval At just prior to t. This is in contrast to Eq. (319) or to the
puff model, for example, which requires an integration over the entire period
t~  t to obtain the concentrations at time t. The second advantage of
Eq. (333) is that it possesses unlimited spatial resolution. As we pointed
out earlier, this property is essential in modeling the microscale concentra
tion distribution. The third advantage of Eq. (333) is that finite differen
cing techniques are not required to evaluate it. In fact, all of the integrals
in this equation can be evaluated analytically when the concentration distribu
tion (c(r,t ,)> is expressed in series form. This advantage circumvents the
> ~ ni /
problems, often associated with differencing techniques, of computational sta
bility and of the imposition of artificial boundary conditions.
2. Development of A Microscale Model of Nonlinear Pollutants Based
on The DiscreteTime, ContinuousSpace Equation
Our aim is to use Eq. (333) as the basis of a microscale model. In this
section, we demonstrate how this might be accomplished, using as the example
a single nonlinear species c whose microscale variations c" are governed by
Eq. (291). The formal relationship between Eqs. (333) and (291) is quite
complicated, and since the details are not important to the present analyses,
we outline the relationship only briefly below.
First, We note that Eq. (219) is of the form
Lc" = kfa + be" + c"2 )+ S" , (334)

169
where
a = c"2 + " , (335a)
b = 2c » (335b)
Consequently, if S" were zero and if c" were initially uniformly distributed
in space (so that LC" = 0), then the chemistry alone would cause c" to change
in time, and we would have
C"(t) = c + q fHr > (336)
LI ~ Q J
where
q = (~c2 _ 4a)1/2 (337a)
»
TTT^J
d  e^LJ3 A ,A J (337b)
C0
and where c!! = c"(0). Recall the form [Eq. (321)] that the kernel p acquires
in situations where the scalar quantity decays at the linear rate
We anticipate, therefore, that if a nonlinear decay of the form found in Eq. (334)
occurs, then the kernel p entering in Eq. (335) can be approximated by
p(r,tr',t') = p'(r,tr',t') [ j + ^ (H~ff)] ' (338)
where, as in the case of Eq. (321), p1 is the solution of the equation
Lp' = 6(r  r') 6(t  t1) . (339)

170
We now state without proof that the equation
dr'
pc(r,tnr',t')
fd +1>
T^Ty
:"(r'
11
dt1 dr1
(340)
becomes equivalent to Eq. (291) in the limit as At = t  t _, approaches zero.
In Eq. (340), q and d are functions of the time and space variables that appear
to the right of the vertical bar in the kernel p. Thus, we can say that Eq.
(340) provides a consistent representation of Eq. (291) inasmuch as the solution
of Eq. (340) can be made arbitrarily close to that of Eq. (291) by making the time
step At suitably small. We should hasten to add, however, that the size of At
is restricted by Eq. (326); so the equivalence of Eqs. (340) and (291) can be
*
realized only by scaling the time axis with some value T and then making At = At/T
suitably small by making T sufficiently large. The implication of this is that
temporal resolution is lost as Eq. (340) is made to approach Eq. (291). Never
theless, for all practical purposes, Eq. (340) will serve well as the working
equation for a model of c".
Consider now the evaluation of Eq. (340). We demonstrate Our proposed
technique for solving this equation by considering only the first integral on
the right side. To simplify the analysis further, we assume homogeneous, sta
tionary turbulence and restrict our attention to a onedimensional problem.
Albeit highly simplified, this case is adequate to demonstrate the essential
features of our technique.

171
Under the conditions that we have imposed upon the turbulence charac
teristics, the kernel p1 becomes
(2*)
'"
r
where
x = x4 + uAt
(341)
(342)
and u and a are functions of x' and t ,.
n1
Turning next to the definition [Eq. (337a) ] of the parameter q that enters
in Eq. (340), we see that.two of the variables upon which it depends, namely,
2
c and c" , spatially averaged quantities. Since the spatial region of con
cern to us in the microscale model is only as large as several grid cells, it
is not unreasonable to treat these two quantities as constants. This assumption
is also supported by observational studies of mean concentration variations in
both the horizontal and vertical directions. The third quantity entering into
r\
the definition of qnamely, <^c' )"is variable over the microscale region,
but for simplicity we also regard it as essentially constant. Combining all of
the above assumptions in Eq. (340), we obtain
c"(x,tn) = 
dx1
where
l
(344a)
a 
(344b)

172
The first term on the right side of Eq. (343) is clearly equal to c(x,t )
[Eq. (333)].
Now let c"(x,t ) be represented by the power series
c"<*. V • ^ an>nx» . (345)
The last term in parentheses in Eq. (343) can now be written in the form
OD
1 "C~"^ m
= j 2^ dnl,mx ' (346)
m  0
where
= Vi,o + c + (347)
(348)
m
d , = aa ,  f > d , , a , . . (349)
n1 ,m n1 ,m 3n /^ n1 ,mk n1 ,k v '
u k = 1
Substituting Eq. (346) into Eq. (343), we obtain
c"(x,tj + c(x,tj = c(x,tj = —iL— , «n_1>m
^0 1^iro m'= 0
• fxmexp( (x X'/ DAt)' }dx'
/
 UAI;
2o*
(350)

173
Changing the integration variable to
5 = x' + uAt  x ,
we obtain
E Vi.
0 u m = 0 •>
(351)
where
x = uAt  x . (352)
Equation (351) can be written in the equivalent form
c(x,t ) = " V™* d , *C~* ryr——ryj I £ 6 d £
m " ° k ~ ° (353)
Evaluating the integral in this expression, we obtain
co
/ * \ M v^ ,1 ^ I/ \m2k 2k
c(x,t ) = £ > d i V m!(x). P ,?. _ ni
1 pn •*«' III,IM ^^ 7 p, \ 17 ?M i *• 'v >
u _ ^iLisy:^iii6.i\y;
m = 0 k = 0
(354)
where
(2k  1)!! = 135.,.(2kl)
and where E1'1^ denotes summation over all k that do not exceed m/2. Upon
expanding the right side of Eq. (354) and collecting like powers of x, and
upon writing the left side of this equation, i.e., c(x,t ), in its series
form
oo
^—^^
f^ [ V "f~ 1 — s f^ V t*f^{Y"f~)
m = 0

174
and grouping these terms with those of the right side of Eq. (354), we obtain
the set of equations for the coefficients d . The desired mircoscale field
c"(x,t ) is now in the form
where
 a
n m = 0 n'm
ansm m
and the f are known functions. The modeling of c" is thus reduced simply to
the algebraic manipulation of coefficients in a series expansion. Mo differ
ential equations or integrals require manipulation.
As we mentioned earlier, this modeling approach has not been developed or
tested heretofore. Consequently, considerable work remains to be done in pro
ducing a working microscale model for photochemical pollutants. We hope,
however, to make considerable progress toward that goal in our continuing model
development program.

175
VI SIMULATION OF BUOYANT PLUMES IN
THE PLANETARY BOUNDARY LAYER
A. INTRODUCTION
Up to this point, all of the work presented in this volume has per
tained to passive scalar quantities, that is, to scalars that do not inter
act dynamically with the fluid. Under this condition, the mathematical
problems of diffusion modeling are greatly simplified because the effects
of the turbulence on the scalar can be described in terms of quantities,
(namely, the turbulent diffusivity) that are properties of the flow. It
is this assumption that made it possible for us in Chapter II to derive
profiles for the diffusivity Kz that were independent of the character
of the diffusing substance.
In air pollution studies involving emissions from power plants, oil
refineries, and other sources in which large quantities of both heat and
pollutants are discharged simultaneously, the "passive scalar" assumption
is not valid on the microscale, i.e., over distances up to about 1 km
from the source. The reason for this is that the heat exhausted along
with the pollutants produces buoyancy forces, which enchance the vertical
transport of the pollutant material. The result is the socalled plume
rise phenomenon, which has been studied analytically by Briggs (1963) and
others. In most of these investigations, attempts were made to describe
the plume rise using semiemm'rical formulas. Some of these formulas
were found to work reasonabl /ell under certain special circumstances,
but not in others. Apparent „ there is no formulation that is applicable
to a wide range of condition^. In this chapter, we outline a model that
we hope will fill this need.

176
The starting point of our simulation is the Lagrangian diffusion
equation introduced in Chapter II [Eq. (16)]:

177
Thus, p is a function of the horizontal projection of the distance to the
release point—but not of the positions of the source or receptor separately
and of the travel time T—but not of the clock time t or release time t'
separately. Consequently, we can calculate p in the following way.
For a given release height z1 , say z1  z,, we release a particle
from each grid point on the horizontal plane at height z, at the initial
instant t when the data begin. If there are N grid points in each hori
zontal plane, then N particles will be released. We allow each particle
to move forward in time under the action of the forces exerted upon it by
the turbulent fluid. By tracking each particle, we create an ensemble of
N particle trajectories, each of duration T, from which the density function
p(£,n,z, z, ,T) , 0 < T < 'T, can be computed by conventional techniques [see,
for example, Bendat and Piersol (1968)]. Here, T denotes the length of
the interval during which data from Deardorff's model are available. In
our early work, T corresponds roughly to 1 hour in real time.
For the function p derived by this process to be applicable
to a given dispersion problem—for example, the diffusion and transport
of pollutants from a stack of height H, diameter D, exhaust velocity V ,
A
and temperature 0 the ensembel of particle trajectories on which the
o
function p is based must adequately reflect the behavior of particles
released under such conditions. In the example just cited, it is clear
that wind shear, buoyancy, and other phenomena have strong effects on
the particle motions. Furthermore, as a result of the proximity of the
particles to one another at their time of release, particle interactions
affect drag forces, heat exchange between the particles, and the atmos
phere, and similar factors. Thus, unless all of these mechanisms are pro
perly modeled in the simulation of the ensemble of particle trajectories,
the density function p based on this ensemble will lead to erroneous es
timates of the mean concentration when used in the <(c> equation.
Hinze (1959) has reveiwed the theory relevant to the prediction of
the motion of discrete, bouyant particles in a turbulent fluid. In the
remaining sections, we use this theory to formulate a model that, using
Deardorff's data as input, can simulate the trajectory ensembeles required

178
to study the diffusion of any type of pollutant (e.g., heat, moisture,
gases) released under any conditions in the planetary boundary layer.
C. MATHEMATICAL BASIS OF THE TRAJECTORY SIMULATION MODEL
PRELIMINARY DESIGN
The Particle Momentum Equation
Let (usv,w) represent the velocity V of the scalar (pollutant) particle
with respect to a fixed Eulerian frame, and let (u ,v ,w ) represent the
average (spatial) velocity V of the fluid in the immediate vicinity of
the given particle, i.e., the local environmental velocity. We distinguish
between the "environmental" velocity V and the ambient fluid velocity Vf,
"** C "" I
because, as a result of the presence of large numbers of scalar particles,
the local velocity of the fluid that a particular particel "sees" will
generally differ from Vf. Also, let T and T represent the particle and
environment temperatures, respectively. Provided that the particles are
sufficiently small, the momentum equation governing their motion can be
written in the form [see Hinze (1959)]
dt
where
[—1
/(a,t)  Ve(x(a,t),t) 
Y Hf Ye(X(a,t)
(356a)
3 =
9v
+
•')
(356b)
(356c)
'(356d)

179
and
r = the particle radius,
p = the particle density,
v = the kinematic viscosity of the environmental fluid,
pg = the density of the environmental fluid,
a = the release points of the particle in question,
X(a,t) = the position of the particle at time t
g  the gravity vector,
and
at • It + v • ' • <357>
The second term on the right side of Eq. (356a) represents a Stokes drag
force on the particle. For this approximation to be valid, the Reynolds
number of the particle with respect to the particle radius r and relative
velocity V  V must be on the order of unity or smaller.
The last term on the right side of Eq. (356a) represents the accelera
tion of the environmental fluid relative to the axis of the moving particle.
In cases such as those of interest to us, where the scalar and fluid par
ticles are dynamically quite similar, accelerations of the one relative to
the other are sufficiently small that the last term in Eq. (356a) can be
neglected. We assume further that the response time a of the particle
is so small, relative to the time scale of changes in V that the quasi
steadystate relation
~nrV(a,t)  0 (358)
holds at all times. Under these two assumptions, we obtain from EQ. (356a)
u = ue , (359a)
v = v , (35Qb)

180
(359c)
In view of the absence of horizontal forces on the particles other than
those exerted by the fluid, we assume that the horizontal components of
the velocities of all particles are equal at all times to the corresponding
velocities of the local ambient fluid. That is, we assume that
u = uf
(360)
v = vf ,
where Uf and Vf are the velocities given by Deardorff's data.
It is more convenient to work with potential temperature than with
densities because the former are given for the fluid by Deardorff's data.
Using the gas law and the definition of potential temperature, i.e.,
R/cp
where p is pressure in millibars, we can convert Eq. (359c) into the
form
(361)
Later, we derive the equation governing the particle temperature 0 that
can be used in conjunction with Eq. (361) to obtain the particle's vertical
velocity component w. However, to render Eq. (361) solvable, we must first
express we in terms of the known ambient fluid velocity component Wf.

181
Toward this end, we consider two limiting cases. First, in the limit
as the travel time becomes large and the initial cloud of scalar particles
becomes so dispersed that each particle is surrounded mainly by ambient
fluid particles, we have
lim wo = wf . (362)
t*» e T
Recall that we represents the velocity of the environmental fluid, com
posed of both scalar and ambient fluid particles, in the immediate vicinity
of the given scalar particle. At the other extreme, where the travel time
T is near zero, the entire collection of scalar particles can be envisaged
as a buoyant jet of diameter D (equal to that of the source diameter)
issuing into the ambient fluid. The velocity of each particle can be as
sumed to be equal to that of the jet and can be approximated by empirical
and theoretical means. For the present, however, we assume that initially
the particles form a spherical cloud of diameter D. This approximation
might apply under very unstable conditions. The initial velocity of euch
particle is now assumed to be equal to that of the entire cloud and is to
be calculated by assuming a drag force on the cloud of the form

3/CD\2
where CD is the drag coefficient of the spherical cloud and uro is the
speed of the cloud relative to the ambient fluid. Invoking the quasi
steadystate assumption as before, we can equate the drag force D and
the buoyancy force B, where
00
\B\ = 9 , (364)

182
to obtain
2
(w  wf)
T
or
w
= wf + sign(s) P^1 . (365)
T JL
Upon substituting this expression into the left side of Eq. (361), we
can find the effective environmental velocity we in the limit as T > 0.
We obtain
4Dg0  0J\1/2 o 2 /0  0
f 
Tim w = w + Sign(e  0   . (366)
This, then, is the environmental velocity component we in the initial
instance where the cloud of scalar particles behaves as a sphere moving
through the ambient fluid. From the two limiting cases given in Eqs. (362)
and (366), we surmise that we can be expressed in the general form
we = wf + <>(T) G(e,0e) , (367)
where
4Dg0  0J\1/2 .2 /9 
G(e,ee) = sign(0  e)  l (368)
and where (0) = 1, Hm (T)  0. (In the general case, G(0,0e) represents
the velocity of the buoyant plume as it emerges from its source.) The rate
at which

183
cloud of particles becomes dispersed,' or, in other terms, how rapidly the
population of particles surrounding any given scalar particle becomes satur
ated with ambient fluid particles. We attempt in a later section to derive
an explicit expression for 4. based on turbulence properties. Thus, from
Eqs. (361) and (367), we arrive at the final form of the expression for w:
? 2
w = wf + <(,(T) G(0,0e) + yL (0 Qe) . (369)
The next task is to derive the equation governing the temperature 0.
2. The Particle Temperature Equation for Dry Plumes
As long as there is a difference between the particle temperature
0 and the ambient fluid temperature 0f, not only is a buoyancy force
exerted on the particle, but also a heat exchange between the particle
and fluid acts to eliminate any temperature differences. This heat flux
can be expressed quantitatively by
q = FA(0  0e) , (370)
where h is the heat transfer coefficient of the particle and A is the
particle's surface area. The temperature change resulting from the heat
flux q is
d0 FA(0  0e)
=~ ' (371)
where Vp is the particle volume and Cp is the specific heat capacity (in
units of energy/mass/°K) . For spheres in the Reynolds number range
1 < Re < 25, F has the form (Kreith, 1966)

184
h  acu"p , (372)
where
^r (2ur)l/
CO
Combining Eq's. (372) and (371), we obtain
(374)
Here it is necessary to express Oe in terms of 0f. For this purpose,
we consider the two limiting cases that we considered earlier in relating
we to Wf. First, in the limit, as the travel time becomes very large, 0e
approaches 0f because each scalar particle "sees" an environment composed
almost totally of ambient fluid particles. We can therefore write
lim ee = 0f . (375)
In the other limit, when T  0, we envisage as before that the entire
collection of scalar particles comprise a spherical cloud of diameter D
and that this cloud exchanges heat with its environment at the same rate
as a rigid sphere would under identical conditions. (In the general case,
we envisage the particles as comprising a buoyant jet, but we do not con
sider this case here.) Moreover, we assume that internal mixing in the
cloud is sufficient to result in an equi partitioning among all of the
particles of the total heat exchanged. In this case, we have
(¥} <«  'f>
limf =h'" ' • 3y , . (376)
pcf'V'"3

185
or
o
dt
(377)
For spheres with Reynolds numbers in the range 25 < Re < 105, the heat
transfer coefficient is given to good approximation by (Kreith, 1966)
h =
0.37kf /u D
~D
0.6
(378)
where kf is the thermal conductivity of the fluid (units of energy/time/
meter/°K). Combining Eqs. (377) and (373), we obtain
1 im
T+0
(0  9f)
dt
1.4 0.6
v
(379)
and by equating this expression and Eq. (374), we obtain the effective
environmental temperature 0e at the initial instant of release:
lim 0 = 0
T+0 e
1 +
2.2rkf(e  0f)
3aPC0(wwf)°4D14v°6
1
(380)
Since both we and 0e are uniquely determined once the populations of
scalar and ambient fluid particles surrounding any given scalar particle
have been prescribed, we assume that the transition of 0e from its initial
value [Eq. (380)] to its final form [Eq. (375)] is describable in terms of
the function cf> used earlier with we. In particular, we assume that
1 +
2.2rkf(0  0f)
1
 (T)]0
f
(381)

186
This equation, together with Eqs. (374), (367), and (369), formsia closed
system of equations that can be solved for the temperature 0 of the par
ticle and vertical velocity component w as a function of travel time T
and the ambient fluid conditions 0f and Wf described by Deardorff's data.
In short, we now have model equations that we can use to simulate the tra
jectories of dry particles of any type in a turbulent fluid. (The deriva
tion of the corresponding equations for vapor plumes is currently in progress.)
All that remains to complete the model is to derive a suitable functional
form for C>(T) .
3. The 4>(i) Equation
As defined, (r) has the limiting values
lim ({.(T) = 1 , (382a)
T+0
lim (T)  0 . (382b)
T+0
In physical terms, <{>(T) represents the fractional population of scalar
particles within a distance of several radii of a given scalar particle.
This suggests that is equivalent to the mean particle concentration in
an expanding cluster. Under the action of turbulence in the inertia! sub
range, a cluster of particles expands at a rate given by
= 20 2/3
30
where e is the energy dissipation rate per mass of fluid, £Q is the ini
tial diameter of the cluster, T is the travel time, and

187
The latter condition may not always be fulfilled in some problems of
interest to us where significant turbulent energy is created by the par
ticle cloud itself. Indeed, observations reveal that under neutral and
stable atmospheric conditions, the turbulent energy within the plumes pro
duced by large power plants is often significantly greater than that in
the atmosphere just outside the plume. Under these_conditions, it will
be necessary to resort to the available empirical and theoretical know
ledge of buoyant jets to derive a saitable approximation for <(>(T). Never
theless, in those cases in which Eq. (383) holds, the mean cloud diameter
d at time T is given by
d(T) . To2 + 10 e2/3 D2/3 rfK _
where D denotes the initial cloud diameter. This expression holds only
as long as d < 10D.
If the initial concentration of particles in the cloud is unity, it
follows from the considerations presented above that
4 D3
I* 8
(385)
We assume that by the time <£/ is outside the inertial subrange, or
d > 10D, and Eq. (385) ceases to hold, either the mignitude of cj> has fallen
to a value near zero or ($2)^ ^ AX, where AX is a measure of the grid size
in Deardorff's model. In the latter case, the growth rate of the cloud can
be estimated explicitly from the computed particle trajectories.

VII SUMMARY
For the reader's convenience, we summarize here the major points
presented in each of the six previous chapters.
A. CHAPTER I
In this chapter, we introduced the term "microscale" to refer to all
phenomena whose space or time scales are too small to be resolvable ex
plicitly and deterministically in grid models of urban pollution. Such
small features result because of turbulence and finite differencing
techniques. Those produced by the former include the turbulent velocity
fluctuations themselves, usually denoted by u1, and the concentration
flucutations c' generated by u1. Both of these microscale features af
fect the spatial and temporal behavior of the evolution of the mean con
centrations of photochemical pollutants. No exact mathematical expres
sions have yet been found that can describe these effects in generalized
situations. Consequently, an important problem in pollution modeling
is the development of an approximate description of these phenomena
that, under most conditions of interest, will be as accurate as the
data upon which the model calculations are to be based. Chapters II,
III, and VI address these problem areas.
The use of finite differencing techniques in diffusion modeling
gives rise to the microscale phenomenon that we have termed the SSV, or
subgridscale concentration variation. This feature occurs because the
spatial variations that exist in the concentration distribution near
point and line sources, such as power plant stacks or highways, are of
much too small a scale to be resolved by the grid meshes used in urban
pollution models. We showed that the SSV can affect the gridaveraged
concentrations of photochemical pollutants in much the same way that the
turbulenceinduced concentration fluctuations can influence the time mean
concentration. We also pointed out that the SSV can complicate the use

189
of pollutant concentration values predicted by a simulation model:
The concentration levels observed at a fixed point will differ from the
spatially averaged concentration predicted at that point as long as the
SSV is not zero. Model validation studies and concentration extrema
forecasts are two examples of applications in which this problem arises.
Chapters IV and V addressed all aspects of the SSV microscale effects.
B. CHAPTER II
In this chapter, we demonstrated how Lagrangian diffusion theory can
be implemented using numerical turbulence models, Using this technique,
we were able to derive expressions for the distribution of the mean con
centration of passive material issuing from a continuous point source
in the planetary boundary layer, We referred to these distributions,
shown in Figures 8(a) and 9(a), as the numericoempirical (NE) solutions
of the Lagrangian equation. Inasmuch as the data sets on which these solu
tions were based were rather small, the results presented in Chapter II
are only tentative. Future studies will attempt to achieve greater ac
curacy and will address the important questions of how the concentra
tions are affected by changes in source height and atmospheric stabilities
(other than those treated here) and by changes to other source types,
such as area or line sources.
Proceeding under the assumption that the accuracies of the present
NE solutions are at worst comparable to those of available empirical
data, we used these solutions as standards for assessing the accuracies
of the three major diffusion models: the diffusion equation and the
Gaussian puff and plume formulas. We devoted subsequent work to optimi
zing the performances of these models by adjusting the functional forms
of the free parameters in each model such that the resulting predictions
were in closest overall agreement with the NE solutions. These "optimal"
parameters, i.e., diffusivities K, and dispersion coefficients a and a ,
Z Z. A
are summarized in functional form in Table 2. With regard to the accura
cies of the optimal models, we note the following:

190
> Of the three models, the diffusion equation is by far the superior
one. Errors are generally on the order of 20 percent, except at
points near the source under neutral conditions in which case much
larger discrepancies are observed [see Figures 8(b) and 9(b)].
> Relative tothe plume formula, the Gaussian puff equation is a
slightly superior model1, but, in quantitative terms, neither
provides an acceptable description of atmospheric diffusion
under neutral stability conditions, at least in the problem
considered in Chapter II. Both models tend to overpredict
groundlevel concentrations arising from elevated sources. Er
rors of 100 percent are prevalent, and in isolated areas they
are much larger. In the particular problem treated, the ac
curacies of the optimal models were acceptable under unstable
conditions (see Figures 18 and 19).
Considering the Gaussian puff and plume models in general, we note
the following points:
> Spatial concentration distributions in the planetary boundary
layer are decidedly nonGaussian.
> The PasquillGifford data commonly used in the plume forumla
are not applicable to emissions from elevated sources unless
some allowance is made for wind shear effects. When such modifi
cations are made, they greatly improve the accuracy of the plume
formula for elevated sources [see Figure 13(b)s as compared with
Figure 18(b)].
The optimal diffusivity, dispersion coefficient, and wind shear pro
files presented in Table 2 were implemented in the form of FORTRAN
function routines. These routines can be used in place of the corres
ponding variable names in existing diffusion models to achieve results
that are compatible with the predictions of Deardorff's boundary layer
model.

191
C. CHAPTER III
Using empirical data, we showed (see Figure 28) that concentration
fluctuations generated by turbulence can dominate the temporal behavior
of the mean concentration of materials that undergo nonlinear chemical
reactions. Although several previous investigators have suggested ap
proximate mathematical expressions for describing these effects, none of
these equations is well suited to pollution modeling studies because
each introduces too many additional differential euqations into the
system to be solved. For this reason, we set for ourselves the task
of developing a new approximation that does not entail multiple equations.
The results of our efforts for the case of a bimolecular reaction are
represented by Eq. (95) for the generalized case, and by Eq. (99) for
the case where the reactants are premixed.
In Section C of Chapter III, we develop for various situations ap
proximate functional forms of the parameters £. and ?R, which enter into
the generalized closure scheme expressed by Eq. (95). [See Eqs. (113),
(121), (124), and (125).] Tests—based on empirical data—of the accuracy
of the overall scheme produced excellent results (see Figures 28 and 30).
D, CHAPTER IV
Just as turbulent concentration fluctuations can affect the mean
rates of nonlinear chemical reactions, so can subgridscale variations
in the concentration fields simulated by numerical models. In Chapter
IV, we derived a test of the significance of these effects [see Eq. (164)],
When this condition is satisfied, SSV effects are negligible and can be
ignored; but, when it is violated, SSV effects may or may not be impor
tant, depending on the particular physical situation. To handle these
cases, we used the concepts that we employed in Chapter III to treat
turbulent concentration fluctuation effects to develop a scheme for
parameterizing the SSV influences. In its most general form, this scheme
is described mathematically by Eq. (182); and in cases where the reac
tants are emitted from the same sources, it takes the form of Eg. (186).

92
We applied the latter to three problems of practical interest:
> A random spatial distribution of point sources, such as building
heating emissions.
> A network of streets of arbitrary separation.
> Strong, isolated point and line sources.
From these applications, we derived a dimensionless number y, defined
by Eq. (232),.that provides a quantitative measure of the magnitude of
SSV effects on nonlinear reactions simulated in urban pollution models.
As portrayed in Figure 34, the magnitude of the SSV impact grows as the
value of y increases. Using representative rate constants, source
strengths, diffusivities, and the like, we evaluated the parameter y for
several of the important photochemical pollutants (see Figure 35).
From this study, we drew the following conclusions:
> The NOOo reaction is the only reaction explicitly treated in the
current SAI model that is affected by subgridscale concentra
tion variations. In this case, the SSV suppresses the effective
rate of ozone depletion.
> Effects from freeways carrying 10 or more vehicles per day are
significant when the wind is parallel to the freeway, but not
when it is perpendicular. Smaller traffic volumes produce pro
portionately smaller effects, and larger volumes produce pro
portionately larger ones.
> SSV effects from networks of city streets carrying 10 or more
vehicles per day are significant when (1) the streets are 150
meters or more apart and (2) the meteorology or source emission
rates are in a state of transition. Under steadystate conditions,
the SSV effects from uniform street networks become negligible,
but those arising from strong isolated sources do not, as described
in the second point above.
These findings are only tentative. In future work, we will attempt
to corroborate them through numerical experiments.

193
E. CHAPTER V
An objection frequently raised against grid point diffusion models
is that in spreading emissions from point and line sources uniformly
throughout the local grid cell, such models produce distorted descrip
tion's of the time concentration field within the immediate vicinity
of point and line sources. This effect complicates the tasks of model
validation and concentration extrema predictions. In Chapter V, we
developed a socalled rnicroscale model that can be used as a subprogram
in any grid model of CO, S02, or other inert or firstorder reactant
to resolve the detailed concentration field at a given point. This sub
model has been implemented in FORTRAN and tested in several preliminary
calculations, as shown in Figure 41.
The chief weaknesses of this microscale model in its present form
are its inability to handle pollutants that react nonlinearly and its
use of Gaussian kernels. The latter precludes applications to problems
of the following types:
> Those in which aerodynamic effects are important, such as
downwash in the lee of an elevated roadway or building.
> Calculations of pollutant levels in street canyons.
> Estimations of pollutant levels on a roadway where vehicle
wake turbulence is primarily responsible for the initial pollu
tant dispersal.
For the purpose of developing a microscale model for pollutants that react
nonlinearly., we derived a special discretetime continuousspace dif
fusion equation [Eq. (333)]. We showed that, when the concentration
field is expressed in a series form and line and point sources are rep
resented in the usual delta function forms, the nonlinear microscale
model reduces to a simple set of algebraic equations (see the last
equation in Chapter V). Tests of this equation and attempts to relax
the Gaussian kernel assumption implicit inthe current microscale model
versions will be undertaken in the next phase of the EPA contract work.

194
F. CHAPTER VI
In this chapter, we outlined a method whereby the boundary layer
model of Deardorff and the equation governing fluid particle dynamics can
be combined to simulate buoyant plumes and cooling tower exhaust in the
planetary boundary later. The idea was to create an ensemble of particle
trajectories from which the probability density function p that enters
into the general Lagrangian diffusion equation [Eq, (16)] can be derived.
The analyses presented in Chapter VI are intended primarily as an illus
tration of the technique we are planning to use to simulate buoyant
plumes. Further theoretical analyses will be required to develop the
final form of the model equations,

195
VIII FUTURE EFFORTS
The work reported in this volume has attempted both to emphasize
the aspects of pollution modeling that are affected by microscale phenomena
and to develop mathematical tools for describing these effects. However,
two basic aspects of this project are still incomplete:
> A thorough and systematic evaluation of the magnitude of
microscale effects in specific problem areas of pollution
modeling interest.
> A thorough testing of the accuracies of the various tools
that we have developed for treating microscale phenomena.
This last deficiency can be remedied by means of a series of straight
forward numerical experiments in which each of the microscale modeling
techniques is examined and exercised under controlled conditions. We
reported the results of a few such studies—the tests of the closure
scheme presented in Chapter III is one example—and others are planned
for the next phase of our work. It is the evaluation of the magnitude
of the actual microscale problem area itself that will require further
planning and careful study. Indeed, the outcome of these evaluations
may well alter the future course of our microscale modeling work. We
foresee three basic steps in the design of these evaluative studies, as
discussed in the sections below.
A. SPECIFICATION OF THE DEGREE OF SPATIAL RESOLUTION
REQUIRED IN URBAN POLLUTION MODELING STUDIES
Air quality standards are currently expressed in terms of certain
timeaveraged concentrations, but no mention is made of the extent of
the spatial areas to which these criteria apply. Considering the fact
that air quality standards are intended primarily as safeguards of public

196
health, one can infer that these standards apply to all points where
plants or animals vulnerable to pollution damage are found. Thus, for
assessments of air quality relative to cropland, livestock, hospital
patients, and so forth, or for validations of a diffusion model using
data gathered by a pollution monitoring station, an urbanscale pollu
tion model possessing point spatial resolution (as described in Chapter
V) should be adequate. However, to assess the pollutant dosages received
by highway patrolmen , bus and truck drivers, road maintenance men.,
and other similar occupational groups, a model should have not only the
capability of providing lineintegrated concentrations, but also the ability
to simulate pollutant concentrations on the sources themselves. The
last feature is outside the scope of present urban diffusion models,
because on the urban scale highways can be treated as line sources of
zero width within which concentrations are infinite. Moreover, on a
highway or within a street canyon, the initial dispersion and the rates
of fast nonlinear chemical reactions are controlled by vehicle wake tur
bulence and other aerodynamic effects that are not described by the con
ventional Gaussian puff and plume models or by the commonly used dif
fusion equation.
These considerations and the likelihood of observing intolerable
highwayintegrated pollutant dosages point to the need for special micro
scale models that can simulate pollutant levels in the motorist's
frame of reference. One of the goals of the next phase of our micro
scale work will be to develop such a model to supplement the capa
bilities of our urbanscale model. We also plan to examine further the
resolution required for various types of modeling studies.
B. EVALUATION OF THE ACTUAL SPATIAL VARIABILITY OF
CONCENTRATIONS AT RECEPTOR SITES OF INTEREST
Having established the locations at which pollutant concentration
calculations are needed and the degree of spatial resolution required at
each site, one is faced next with the following question: Are local

197
spatial variations in the concentration field so large that an urban
scale model alone is incapable of providing a representative estimate
of the true pollutant levels at those sites? In the case of the fixed
point dosage estimates, this question is essentially one of whether local
sources are responsible for a significant fraction of the pollution ob
served; and in the case of the roadwayintegrated dosage calculations,
the question is whether background concentrations are so smallcompared
with those on and produced by the roadway—that the urbanscale model
itself is needed.
Continuous pollutant measurement data can help resolve these ques
tions. For example, suppose that such a record is available for a
groundlevel traverse of a city. If moving averages of this record were
made to remove all variations that a grid model of that pollutant could
resolve, then the residual concentrations would represent the micro
scale variations in question. If these microscale deviations in populous
areas turn out to be only a small fraction of the total concentration
observed at those same locations, then urbanscale models alone, i.e.,
without a supplementary microscale module, should be adequate for assessing
urban air quality.
Analyses of the type just described have been performed on unpub
lished oxidant measurements made by Lamb and Neiburger in the Los
Angeles basin in 1967. Preliminary results of this work indicate that
significant microscale oxidant variations occur on freeways themselves.
However, the amplitude of these perturbations decreases very rapidly with
distance, both upwind and downwind, from the freeway. Microscale varia
tions in oxidant levels also occur on heavily traveled city streets and
in street canyons, but the amplitudes are less than those observed on
freeways. Finally, on all streets with little or no traffic, microscale
oxidant variations are negligible.

These preliminary findings emphasize the need for a microscale model
applicable to roadways themselves. They also suggest that, except for
sites on streets or freeways, it should be possible to obtain accurate
fixedpoint assessments of oxidant levels using an urbanscale model
alone (one in which subgridscale concentration variations have been prop
erly parameterized). We plan to continue these empirical studies so
that we can determine where the true microscale modeling problems lie.
C. ASSESSMENT OF THE NEED FOR REFINED MICROSCALE
TRANSPORT AND DIFFUSION FORMULAS
The foregoing discussions point toward the need for pollution simu
lation models applicable to both freeway and street canyon environments.
Since within such regions transport, diffusion, and concentration fluc
tuation effects are controlled by vehicle wake turbulence and other aero
dynamic effects not normally considered in conventional pollution modeling
theories, the question arises of whether refined diffusion formulations
will be required to develop the microscale models needed. We have al
ready begun to develop a model of wake turbulence that can be used in
conjunction with the theory presented in Chapter III of this report to
simulate chemical reactions in the wakes of motor vehicles. We also
plan to develop expressions for the kernel p, which enters in the
Lagrangian diffusion equation [see, for example, Eq. (295)], to render
this equation applicable to situations where aerodynamic effects, such
as building and elevated roadway wakes, are important. This work will
be described in subsequent phases of our EPA contract effort.

199
APPENDIX A
FORTRAN PROGRAMS FOR USTAR, DKZ, AND UBAR

USTAR =
ROUTINE
U10 = wind speed at anemometer height,
Z10 = anemometer height, in meters,
ZI = inversion height, in meters,
ZO = surface roughness, in meters,
ZI0VL = ZI/L = stability parameter.
FUNCTION USTAR (U10,ZO,ZI0VL,ZI,Z10)
IF (ZI0VL) 10,20,30
10 Z10=Z10/ZI
ZO=ZO/ZI
IF (Z10.GT.0.025.AND.ZO,LT.0.004) GO TO 15
X=(1.15.*Z10+ZO)*ZI0VL)**0.25
XO=(1.15.*ZO*ZI0VL)**0.25
A1=ATAN(X)ATAN(XO)
A2=ALOG((X1.)/(XO1.))ALOG((X+1.)/(XO+1.))
U10BAR=(2.*AHA2)/0.35
USTAR=U10/U10BAR
RETURN
15
RETURN
20
IF (Z10.GT.0.3)Z10=0.3
UU26.22+153. 2*Z1014.28. *Z10**2+5541. *Z10**3
7523.*Z10**4ALOG(ZO*6.8E6)/0.35
USTARU10/UU
USTAR=0.35*U10/ALOG((Z10+ZO)/ZO)
RETURN
C EXPRESSIONS BELOW FROM RAGLAND PAPER
30 ZZ=Z10*ZI0VL/ZI
IF (ZZ.GT.1.0) GO TO 35
USTAR=0.35*U10/(ALOG((Z10+ZO)/ZO)+5.2*ZZ)
RETURN
35 USTAR=0.35*U10/(ALOG((Z10+ZO/ZO)+5.2)
RETURN
END

VERTICAL DIFFUSIVITY ROUTINE
Z = height of level where KZ is wanted,
ZI = inversion height,
USTAR = friction velocity,
F = Coriolis parameter = 2ft sin 0,
ZI0VL = ZI/L = stability parameter.
FUNCTION DKZ(Z,ZI,USTAR,F.ZI0VL)
IF (ZI0VL) 10,20,30
10
RETURN
15
RETURN
20
' RETURN
25
RETURN
30
35
RETURN
RETURN
END
Z=Z/ZI
IF (Z.GT.1.0) GO TO 15
DK=6.934E3+0.6H3*Z+3.297*Z**2
6.442*Z**3+3.153*Z**4
IF (DK.LT.0.0) DK=0.0
DKZ=DK*USTAR*ZI
DKZ=0.6123*EXP((Z1.)**2/.02)
Z=Z*F/USTAR
IF (Z.GT.0.45) GO TO 25
DK=7.396E4+6.082E2*Z+2.532*Z*Z
12.72*Z**3+15.17*Z**4
DKZ=DK*USTAR*USTAR/F
DK=3.793E3*EXP((Z0.45)**2/2./4.E2)
DKZ=DK*USTAR*USTAR/F
XL=ZI/ZIOVL
IF (Z.GT.085*XL) GO TO 35
DKZ=.35*USTAR*Z/(1 .+4.7*Z/XL)
DKZ=0.0

202
WIND SPEED PROFILE ROUTINE
FUNCTION UBAR (Z,ZO,ZI,ZI0VL,F,USTAR)
IF (ZI0VL) 10,20,30
10 ZO=ZO/ZI
2=1/21
IF (Z,GT.0.025.AND.ZO.GT.0.004) GO TO 14
X=(1.15.*(Z+ZO)*ZI0VL)**0.25
XO=(1.15.*ZO*ZI0VL)**0.25
A1=ATAN(X)ATAN(XO)
A2=ALOG((X1.)/(XO1.))ALOG((X+1.)/(XO+1.))
UBAR=(2.*AHA2)*USTAR/0.4
RETURN
14
RETURN
15
RETURN
20
RETURN
25
RETURN
30
RETURN
35
IF (Z.GT.0.3) GO TO 15
UU=26.22+153.2*Z1428.*Z**2+5541.*Z**3
7523.*Z**4ALOG(Z0*6.8E6)/0.35
UBAR=UU*USTAR
UBAR=USTAR*(32.33ALOG(ZO*6..8E6)/0.35)
ZO=ZO*F/USTAR
Z=Z*F/USTAR
IF (Z.GT.0.055.AND.ZO.LT.0.006) TO TO 25
UBAR=USTAR*ALOG((Z+ZO)/ZO)/0.37
IF (Z.GT.0.21) Z=0.21
UU29.82+213.2*Z1989.*Z**2+8743.*Z**3
14670.*Z**4ALOG(ZO*1.5E7)/0.37
UBAR=UU*USTAR
ZZ=Z*ZI0VL/ZI
IF (ZZ.GT.1.0) GO TO 35
UBAR=USTAR*(ALOG((Z+ZO)/ZO)+5.2*ZZ)/0.37
IF (ZZ.GT.6.0) ZZ=6.0
Z=ZZ*ZI/ZI0VL
WRITE (6,1000)
UBAR=USTAR*(ALOG((Z+ZO)/ZO)+5.2)/0.37
RETURN
1000 FORMAT (1HO, 'NO ACCURATE WIND DATA ABOVE SURFACE LAYER IN STABLE CASE')
END

203
APPENDIX B
DERIVATION OF EQS, (193) AND (195)

204
APPENDIX B
DERIVATION OF EQS, (193) AND (195)
Equation (193) states that
Aj(r,t) = C f p(r  r',t,t') S(r',t') dt1 dr'
0
and Eq. (195) reads as follows:
A?(r,t) =fff f p(rr',t,t') p(rr",t,t") G(r' ,r" ,f ,t") dt1 dt" dr1 dr"
1 ~ I / / I ~ ~ ~~ ~~
JJ JQ JQ
In the derivation of these equations, we first want to show that
l) d' d =
'J
Av(r)
(Bl)
where
" xAx Ay AZ
Ax Ay AZ
(B2)
Since each of the three integrals in /dr' in Eq. (Bl) are similar, we con
sider the onedimensional case only, i.e.,
/X+AX r X+AX •
J J P(x"  x1) C(x") D(x') dx1 dx" = / / p(x"  x1) C(x") D(X') dx" dx1 .

205
Let
x" = x  x1 + x1 =>dx" = dx1
Then
X TAX
f f p(x  xj) C(x  x] + x1) D(x') dx] dx
w ' A v
y.X.+AX
/ C(x  x] + x') D(x') dx1 dx
]
X,AX
Now let
t. = x'  x, $> dx1 =
r r Ax
= / P(x  x^ / C(x + 0 D(x] + 0 d? dx]
Ax
AX
Upon repeating for the y and z integrals, we obtain Eq. (Bl)
Now we want to show that
1 f \ ff p(r  r1) p(r  r") S(r') S(r") dr1 dr" dr
v J \JJ ~ ~ ~~ ~  ~ ^ j 
AV L
= f f P(r  ^') P(r  r") S(r') S(r") dr1 dr"
^yy ~~
where the product mean value in the right integral is defined as in Eq. (B2)
To prove this, we start with the onedimensional problem as before:

206
x. X+&X
J ff
 x'} p(c  x"} S(x') S(x") dx' dx1'
XAX
l rx+*x r
= 2& J y P(C  x') S(x') F(c)
dx'
(B3)
XAX
where F(c) = p(c  x")S(x")dx". We can use Eq. (Bl) to write Eq. (B3)
in the form
P(C  x1) S(x') F(s) dx1
(B4)
where
S(x')
Now let
rAx /•
J JpU  x" + 0 S(xr + 5) S(x") dx"
AX
x" = X" + r => dx"  dx"
Then
S(x') F(c) =
2A
 x") S(x' + 0 S(x" + 0 dx"
AX
= f P(C  x") S(x') S(x")
dx"
where the tilde average here is defined as in Eq. (B2). Using this result
in Eq. (B4), we obtain

207
/ iff P(~ " ~'} p(r  r"} S

208
APPENDIX C
FORTRAN LISTINGS OF THE CALC, SIGMAX,
AND SIGMAZ FUNCTIONS

209
SUUROUTINE CALC(XALPHA,YALPHA,ZALPHA.XR„YR,2R.ISTYPE.XLNGTH,THETA»
V DELTAX.D EL TAY,DELTAZ,UB/iR,VeAR,RESULT)
COMMON/M ICF<0/ IT F MAX , NCHFCK. N S A MR , X S AMP , F PSLON , I X , R I TF . S I GF A C
co MM ON/T A RLE:/ ERFTC isos) ,EXPT{2005 ) .NERF,NEXP,ERFLIM,DERF.DEXP,
^ . tz XPL I M — • •   ..._T  _,T _ . ._
DIMENSION TEST(SO)
LOGICAL HALT,RITE.CALL1
IOUT = ITRMAX/NCHECK •«• 1
SQ2 =—1 ,4 14 21 35 : —
SQ2PI = 2.5066283
T2PI = 6.2821653
DELTA0.0
.___ IF_{ I STYPE ) 20,10,15
10 COSTH = COS(THETTA)
SINTH = SIN(THETA)
XLO2 = XLNGTH/2.
DELTA=XLC2
GO TO 20
15 DXO2 = OELTAX/2
DYO2 = DELTAY/2
2= DELTAZ/2,
OELTVe  DELTAX*DELT AYfOELT AZ*8
DELTA=DXO2
20 XMXA = XPXALPHA
_____ I _ yMYA_=_YRY ALPHA ___
ZMZA = ZRZ^LPHA
' ZPZA = ZRtZALPHA
ZMZA22 = ZMZA«*2/2.
ZPZA22 = ZPZA**2/2.
R = SORT (XMXA<XMXA+YMYA*YMYA)
VEL = SQRTC UBARwOBARf VBAR*VBAR)
TO = R/VEL
SIGX=SJGf/AX(TO)
SGDELTSIGFAC*SIGX+DELTA
TUP=(R+SGDELT)/VCL
TLO=

210
60 EXPARG=(XUHART*XU0ART»YVBART*YVBART)/2.0/SIGX?
PXPY=EXPGCEXPARG)/SIGX2/T2PI
65
GO TO 75
70 Al = < XUBART+DX02) /SQ2SGX
• A2=(XUBARTDXC2 >/S02SGX
= ( YV8ARTSDYC2) /S02SGX _
^= (YVBAKTDY02 )/SQ2SGX
* PXPY= = T22
IF( ITROUTNSAMP.LT.O) GO TO 200
X  CONVERGENCE ROUT I NE.
SUMJ=0.0
JMIN=1TROUTNSAMPM
DO 150 J=JM! K, ITRCUT
 150  SUMJ = SUM J«TEST(J ) _____
SBAR=SUV J/XSAMP
SUMJ=0. 0
DO 160 J=JMIN, ITROUT
.4..60  SUMJ = SUMJ+ { TFST( J) SBAR )**2.
HALT=SQRT C SUMJ/XSAMP).LE.EPSLON*SOAR
IF(HALT) GC TO 205
CONTINUE
= 7RNGE*T22
IF{.NOT.RITE) PETURN
DO 280 JTEST=1,ITROUT
ITRNJTEST^NCHECK
280 WRITE (6, 1042) —ITRN,TESTfJTEST)
RETURN
1028 FORMAT(1H0.7E13.5 )
1030 FORMAT(1H00'RECEPTOR LIES ON A POINT SOURCE')
1042 FORM AT UH.J AFTER' » I6.2X. MTERAT IONSTEST= •» El 2. 5)
END
SUBROUTINE TABLET
— COMMON/TABLE/£RFTt 1 505 ),EXPT( 2O05 ) .NERF.NEXPsERFUIH^DERF^OE XP,
* EXPLIM
Y=0.0
N02=NERF/2
DERF=ERFLI M/FLOAT(N02)
ERFTCN02P1 )=0 .0
DO 15 N=1,NO2
ERE=ERF( Y)
• ERFT( NQ2P1+N) =ERE
IS ERFT{N02P1N) =ERE
   DEXP=EXPLJ
X=DEXP
RETURN
END
DO 20 N=1,NEXP
X = X + DE:XP
EXPT=EXP<X)

211
FUNCTION EXPO(APG)
^COMMON/TABLE/ ERFT(1505),£XPT<2005KNERF..NEXP.ERFLIM.OERF^OEXP.
IE=IFIX(/»RG/DEXP)+1
IF(I5.GT.NFXP) IE=NEXP
 EXPO=EXPT< IE) 
RETURN 
END
^FUNCTION ERR OR (ARC) .
^COMMON/TABLE/ ERFT ( 1 505 ) , EXPT ( 20 05 ) . NERF , NEXP
RETURN
END—
EXPL1M
IERF=IF IX( (ARG+EPFL IM)/DERF
IF( IERF.LT»t).IERF=1
IF( IERF.GT .NERF) IERF = NERF
ERROR=ERFT(IERF)
TN 5 T_GM A7. ( T~P
COMMON WS.TAR,USTAR,H,F
10 TS= T*WSTAR/H
IF(TS.GT.Oc525) GO TO 15
SIGMAZ=H*( 0 .65*TS1 . 89*TS*TS +5
PFTIIPM ,
15 IFCTS.GT.l .1 ) GO TO 16
SIGMAZ=H*(1. 11 1*TS*TS + 2.344*TS0.33.4
RETURN
. RETURN
20 TS=T*F/0.45
USTARF=LSTAR/F
 LF(TS.GT^O.OS) GO TG2 ______
SIGMAZ=0.45*USTARF*0.64*TS
RETURN
21 IF(TS.GT.0.175) GO TO 22
SIGMA.Z=X1 *45JJUSTJXRF*< 0 .Ol
RETURN
22 SIGMAZ=0.45*USTAPF*( 0.058*0 .12"9*TS)
IF( SIGMAZ. GT.0.31*USTARF) S I GM AZ=0 . 3 1 *UST APF
LIMIT nnr^ rn M A x ( *.\r,^/u ) = .Q IN IJM^TAPI F r.A.sF AND H= .
NEUTRAL CASE.
RETURN
30 WRITE (6, 1003)
RETURN
1003 FORMATC 1HO.. 'NO SIGMA DATA FOR STABLE CASE1 >
END
TI 0 N S I G M A X ( T f t
COMMON"~wsT~A~Dt USTAK.H. F
IF(WSTAR) 30,20.10
__TS=TJf WSTAR/H
IF(TS.GT.O.OS) GO TO 15
SIGMAX=2.2*TS*H
RETURN
15
RETURN
.20 TS=T*F/0.45
IF(TS.GT.O .1 ) GO TO 22
RETURN
22 SIGMAX=( .45*USTAR/F)*(0.1»0.92*TS)
RETURN
 30  :  WRLTE (6. 1 003 X 
S1GMAX=0.5*(USTAR/F)
RETURN
1003 FORMATC 1HO» "NO DATA FOR STABLE CASE*)
_ END _______________ .

212
APPENDIX D
THE MONTE CARLO TECHNIQUE
USED BY THE CALC SUBPROGRAM

213
APPENDIX D
THE MONTE CARLO TECHNIQUE
USED BY THE CALC SUBPROGRAM
We propose to use a Monte Carlo technique to evaluate the integrals
entering in Eqs. (302), (305), and (306). The technique is best des
cribed by the mathematical analyses that are required to prove its validity.
Consider a function f(t) defined in the interval t« < t 5 t + T
as shown in Figure Dl. Suppose we pick a sequence of numbers t.,
i = l,2,...,n, at random in the interval tQ < t 1 tQ + T. We choose the
numbers by a random process such that each number in the interval is
equally likely to occur. For eacl
a corresponding unique number f.:
equally likely to occur. For each number t. in the sequence, there is
f. = f(t.O
i v i
Thus, the probability of observing a value f in the range f. < f  f. + f
is by definition
m(e.)
p(f.) = prob(fi < f < f. + f) = ^ , (Dl)
where m(e.) is the measure of the set e. in which f. < f(t)  f. + f
(see Figure D2) in the domain (tQ, tg+T).

f(t)
t0+T
FIGURE D1. A FUNCTION f(t) DEFINED IN THE INTERVAL tQ < t <_ tQ + T

f(t)
A
fi
TOTAL LENGTH = m(e.) = MEASURE OF SET e.
FIGURE D2. ILLUSTRATION OF THE PROBABILITY THAT f(t) OCCURS
IN THE INTERVAL fn TO fj + Af OVER THE TIME t0
VT
rv>
en

21G
Now, by definition of the Lebesque integral, we have
,t0+T
f(t) dt 
lim
Af>0
(D2)
where the summation is over all intervals f. in the range °° < f < °°.
J
Let the mean value of the random sequence f(tj), i = 1,2,... ,m, generated
above be noted by
(D3)
Then by definition the ensemble mean value of the random variable f(tj) is
f = lim f = lim 2,
" Af+0 j
(D4)
where p(fj) = prob(fj < f <_ f j + Af) and the summation is over all
intervals fj in the infinite range of f values. From Eqs. (Dl) and (317),
we obtain
f = lim 2, f;
Af+0 j J
m(e.
T
Comparing this with Eq. (D2), we see finally that
f
,+T
f(t) dt = T li
m
n
1 2, f(ti
(D5)
on
This result means that we can approximate the integral of any functit
f(t) over any domain by simply multiplying the length of the domain T by
the mean value of the random variable f(tj) formed by picking points tj
at random in the domain of integration.

217
The speed of execution of this Monte Carlo integration technique
can be increased greatly by tabulating the exponential and error func
tions that appear in the integrands of several of the integrals of
interest. That is, rather than use the EXP and ERF FORTRAN routines to
compute the function value each time it is required, vie create a table
of the values of each function initially and look the value up in the
table as needed. This procedure has been found to reduce computation
time by about onethird. The subroutines and functions TABLET, EXPO,
and ERROR listed in the CALC program replace the EXP and ERF functions.

218
APPENDIX E
)F rAB, rA,
MULTIJET REACTOR DATA AND TOOR'S THEORY
DERIVATION OF PAB, i?A, AND ffi FROM

219
APPENDIX E
DERIVATION OF TAB, JT, AND TR FROM
MULTIJET REACTOR DATA AND TOOR'S THEORY
The purpose of this appendix is to give a brief description of the
multijet reactor data and how we used the data to calculate the inert
/\ <*. A
correlation functions (r*n, l\, and rR) used in Chapter III.
As noted in Chapter III, Toor (1962) developed a theory that relates
the concentration statistics of two species undergoing a very rapid irre
versible reaction in a turbulent mixer to the concentrations of inert
species in an identical mixer. To test the theory experimentally, Vassilatos
and Toor (1965) designed an ideal onedimensional tubular reactor having a
head made of some 100 small nozzles (Figure El). Reactants are fed through
alternate jets to simulate a crosssectionally uniform concentration profile.
A modified version of this reactor was later made and used by Mao (1969).
Extensive measurements covering a wide range of reaction speeds were made
in these reactors. However, the results provided only a partial test of
Toor's theory because no inert mixing data were taken. Later measurements
of inert concentration statistics in the same reactors by McKelvey (1968)
and Zakanycz (1971) have corroborated Toor's theory.
polymethyl methacrylate
epox.y
:1.5
4.5
FIGURE El. SIDE AND END VIEWS OF TUBULAR REACTOR

220
The theory presented by Toor (1962) states that for very rapid reactions
with stoichiometric feed,
*
, , 2
where x = t is the axial position along the reactor length and XQ is
some reference point where the reactor has reached the state of cross
sectional homogeneity. Toor and his coworkers designed the multijet reactor
so that homogeneity is achieved virtually at the inlet of the reactor. Thus,
we set x~ = 0. Toor (1969) further derived a relationship between the con
centration fluctuations of two unpremixed inert species fed into this reactor:
(E2)
2
(E3)
Since the species are not premixed,
=  . (E4)
Combining Eqs. (El) through (E4), we obtain