-------
«»•
^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) d-c 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(r-r',t,t') p(r-r'-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 subgrid-scale
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-
of-magnitude estimate of the effect of the subgrid-scale concentration variations
on the chemistry.
According to our scheme, the measure of the importance of the subgrid-scale
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 subgrid-scale 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
Subgrid-scale effects will be most pronounced for fast reactions in which T
7
is small. For the case of the NO-0- 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, subgrid-scale 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 subgrid-scale 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 i-c • — / I l^.£_E * *
where Q is the emission rate,
2 ,, x
°z - kz - '
a2 = K * ,
y y u
and where the x-axis lies on the plume center! ine with the mean wind parallel
to the x-axis. 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 subgrid-scale 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 NO-0- 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 subgrid-scale effect diminishes.
These results agree with the analyses of randomly distributed point sources
presented in Section IV-E-1.
Turning finally to the subgrid-scale 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 subgrid-scale-induced perturbation
in the concentration will be on the order of
-1/2
u AZ
c 1-AB 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 NO-CL 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 subgrid-scale chemistry
effect occurs.
We conclude that in contrast to the point sources, line sources of commonly
encountered strengths produce negligible subgrid-scale chemistry effects, at
least insofar as the N0-03 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 subgrid-scale 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
urban-scale 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 small-scale 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 small-scale 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 urban-scale 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 urban-scale model. That is,
we sought the capability of examining the near-source small-scale concentra-
tion distribution at any point without having to make multiple runs of the
large-scale model. We also wanted the microscale model not to be structured
around or coupled to the urban-scale 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 i-th 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 error-free 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 small-scale
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)
X-AX y-Ay Z-AZ
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 subgrid-scale 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 subgrid-scale variations
c" can arise only from subgrid-scale 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 subgrid-scale field c". Thus, our approach
to the restoration of point resolution to the grid model is simply to develop
a microscale, or subgrid-scale, 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 subgrid-scale 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 SUBGRID-SCALE VARIATIONS c"(r,t)
For simplicity, we consider a one-dimensional problem governed by the
equation
2
|f = K-^-|+ fi(z) 6(t) (2-75)
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 f-z
4Kt
(277)
Now let
z+Az
c(z,t) =
2Az
c(z',t) dz1
Z-AZ
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
subgrid-scale 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 SUBGRID-SCALE CONCENTRATION VARIATIONS
ARISING FROM AN INSTANTANEOUS POINT SOURCE
-------
145
That is, o* is the approximate half-width 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 steady-state 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 x-axis). Translated
into terms of the two-dimensional 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 second-order 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 subgrid-scale 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 subgrid-scale 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 subgrid-and supergrid-scale 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 subgrid-scale 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,t|r't') exp |~-k(t - t1)] S"(r',t') dt1 dr1
'O
+ f p(r,t|r'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 a-th point source (of
pa
which N are present), r is the position of the a-th 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 a-th point source, L. is the con-
a p
tribution to c" from the 3-th line source, and where V is the collective con-
Y
tribution of all sources in the y-th grid cell. Below we derive the functional
form of each of these terms, assuming that in the open terrain and quasi-steady
flow regimes of interest to us here the kernel p in Eq. (295) has the Gaussian
form
-------
150
p(r,t|r',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 low-level 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 3-th 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 y-axis 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 urban-scale 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)
Z-j _____ 7-4-7 -4-
-1 - O i. ~ i. ~ ~7T
B =
and (x ,y , z ) denotes the center of the y-th 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 I-th point source.
ZPS(I))
XLS(J) ) _ (x,y) coordinates of the center of the J-th line
YLS(J) | line source (ZLS = 0 assumed).
THETA(J) = angle of the J-th line source (see Figure 38).
XLNGTH(J) = length of the J-th line source.
SP(I) = strength (mass/time) of the I-th point source.
SL(J) = strength (mass/length/time) of the J-th 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 self-explanatory, 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 full-scale 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
10-METER
STACK
.GROUND-LEVEL
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
point-scale 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 V-B 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 wedge-shaped 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 V-B 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 wedge-shaped
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 A-A 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 i-th 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 line-source 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 i-th
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 print-out 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 SUBGRID-SCALE
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 o-f 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 subgrid-scale problem:
because due to the presence in S" of delta functions arising from point and
line sources, subgrid-scale 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
well-known 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 Discrete-Time, Continuous-Space 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,t|r',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,t|r'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 discrete-time,
continuous-space 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'~n-l ~n-2' '~0 K ~n ~n-l ' 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 Chapman-Kolmogoroff 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 so-called
energy-containing 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 Chapman-Kolmogoroff 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
/-/n-l
' 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,
\ ~ n-1 '
Eq. (327) is equivalent to
-------
168
^p(r,tn|r',tn_1)dr'
P(r,tn r',t') S(r',t') dt1 dr' . (333)
'n-1
This is the discrete-time, continuous-space 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 t-At 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
> ~ n-i /
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 Discrete-Time, Continuous-Space 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 f-Hr > (336)
LI ~ Q J
where
q = (~c2 _ 4a)1/2 (337a)
»
TTT^-J
d - e-^L|J3 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,t|r',t') = p'(r,t|r',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,tn|r',t')
fd +1>
T^Ty
:"(r'
1-1
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 one-dimensional 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 ,.
n-1
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 q--namely, <^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^ dn-l,mx ' (346)
m - 0
where
= Vi,o + c + (347)
(348)
m
d , = aa , - f- > d , , a , . . (349)
n-1 ,m n-1 ,m 3n /^ n-1 ,m-k n-1 ,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
• fx-mexp(- (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-—-—r-yj I £ 6 d £
m " ° k ~ ° (353)
Evaluating the integral in this expression, we obtain
co
/ * \ M v^ ,1 ^ I/ -\m-2k 2k
c(x,t ) = -£- > d i V m!(-x). P ,?. _ ni|
1 pn •*-«' II-I,IM ^^ 7 p, \ 17 ?M i *• 'v-- >
u _ ^iLisy:^iii6.i\y;
m = 0 k = 0
(354)
where
(2k - 1)!! = 1-3-5.,.(2k-l)
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 so-called 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 semi-emm'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-
steady-state 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-
steady-state assumption as before, we can equate the drag force D and
the buoyancy force B, where
0-0
\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
4Dg|0 - 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
4Dg|0 - 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(w-wf)°-4D1-4v°-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
subgrid-scale 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 grid-averaged
concentrations of photochemical pollutants in much the same way that the
turbulence-induced 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 numerico-empirical (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 to-the 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
ground-level 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 non-Gaussian.
> The Pasquill-Gifford 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 subgrid-scale 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 NO-Oo reaction is the only reaction explicitly treated in the
current SAI model that is affected by subgrid-scale 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 steady-state 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 so-called rnicroscale model that can be used as a subprogram
in any grid model of CO, S02, or other inert or first-order 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 discrete-time continuous-space 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 in-the 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
time-averaged 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 urban-scale 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 line-integrated 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
highway-integrated 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 urban-scale 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 roadway-integrated dosage calculations,
the question is whether background concentrations are so small--compared
with those on and produced by the roadway—that the urban-scale 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
ground-level 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 urban-scale 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
fixed-point assessments of oxidant levels using an urban-scale model
alone (one in which subgrid-scale 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((X-1.)/(XO-1.))-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
UU-26.22+153. 2*Z10-14.28. *Z10**2+5541. *Z10**3
-7523.*Z10**4-ALOG(ZO*6.8E6)/0.35
USTAR-U10/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.934E-3+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(-(Z-1.)**2/.02)
Z=Z*F/USTAR
IF (Z.GT.0.45) GO TO 25
DK=7.396E-4+6.082E-2*Z+2.532*Z*Z
-12.72*Z**3+15.17*Z**4
DKZ=DK*USTAR*USTAR/F
DK=3.793E-3*EXP(-(Z-0.45)**2/2./4.E-2)
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((X-1.)/(XO-1.))-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*Z-1428.*Z**2+5541.*Z**3
-7523.*Z**4-ALOG(Z0*6.8E6)/0.35
UBAR=UU*USTAR
UBAR=USTAR*(32.33-ALOG(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
UU-29.82+213.2*Z-1989.*Z**2+8743.*Z**3
-14670.*Z**4-ALOG(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(r-r',t,t') p(r-r",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)
(B-l)
where
" x-Ax Ay AZ
-Ax -Ay -AZ
(B-2)
Since each of the three integrals in /dr' in Eq. (B-l) are similar, we con-
sider the one-dimensional 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 - x-j) 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. (B-l)
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. -(B-2)
To prove this, we start with the one-dimensional problem as before:
-------
206
x. X+&X
J ff
- x'} p(c - x"} S(x') S(x") dx' dx1'
X-AX
l rx+*x r
= 2& J y P(C - x') S(x') F(c)
dx'
(B-3)
X-AX
where F(c) = p(c - x")S(x")dx". We can use Eq. (B-l) to write Eq. (B-3)
in the form
P(C - x1) S(x') F(s) dx1
(B-4)
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. (B-2). Using this result
in Eq. (B-4), 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 3-5- : —
SQ2PI = 2.5066283
T2PI = 6.2821653
DELTA-0.0
.___ IF_{ I STYPE )- 20,10,15-
10 COSTH = COS(THETTA)
SINTH = SIN(THETA)
XLO2 = XLNGTH/2.
-DELT-A=XLC2
GO TO 20
15 DXO2 = OELTAX/2
DYO2 = DELTAY/2
2-= DELTAZ/2,
OELTVe - DELTAX*DELT AYfOELT AZ*8
DELTA=DXO2
20 XMXA = XP-XALPHA
_____ I _ yMYA_=_YR-Y ALPHA ___
ZMZA = ZR-Z^LPHA
' ZPZA = ZR-t-ZALPHA
ZMZA22 = ZMZA«*2/2.
ZPZA22 = ZPZA**2/2.
R = SORT (XMXA<--XMXA+YMYA*YMYA)
VEL = SQRTC UBARwOBAR-f VBAR*VBAR)
TO = R/VEL
SIGX=SJGf/AX(TO)
SGDELT-SIGFAC*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=(XUBART-DXC2 >/S02SGX
= ( YV8ART-S-DYC2) /S02SGX _
^= (YVBAKT-DY02 )/SQ2SGX
* PXPY= = T22
IF( ITROUT-NSAMP.LT.O) GO TO 200
X -- CONVERGENCE- ROUT I NE.
SUMJ=0.0
JMIN=1TROUT-NSAMP-M
DO 150 J=JM! K, ITRCUT
-- 150 -- SUM-J = 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
ITRN-JTEST^NCHECK
280 WR-I-TE (-6, 1-04-2-) —I-TRN,T-EST-fJTES-T-)
RETURN
1028 FORMAT(1H0.7E13.5 )
1030 FORMAT(1H00'RECEPTOR LIES ON A POINT SOURCE')
104-2- FORM AT UH--.J AFTER' » I6.2X-. MTERAT IONS-TEST= •-» El 2. 5)
END
SUBROUTINE TABLET
— COMMON/TABLE/-£RFTt 1 505 )-,EXP-T-( 2O05 ) .NERF-.NEXP-sERFUI-H^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{N02P1-N) =-ERE
- - - DEXP=EXPL-J
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*TS-1 . 89*TS*TS +5
PFTIIPM ,
15 IFCTS.GT.l .1 ) GO TO 16
SIGMAZ=H*(-1. 11 1*TS*TS + 2.344*TS-0.33.4
RETURN
. RETURN
20 TS=T*F/0.45
USTARF=LSTAR/F
- LF(TS.GT^O.OS) -GO -TG-2 ______
SIGMAZ=0.45*USTARF*0.64*TS
RETURN
21 IF(TS.GT.0.175) GO TO 22
SIGMA.Z=X1 *45JJ-USTJXRF*< 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 - : -- WR-LTE (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 D-l. 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) = -^ , (D-l)
where m(e.) is the measure of the set e. in which f. < f(t) - f. + f
(see Figure D-2) in the domain (tQ, tg+T).
-------
f(t)
t0+T
FIGURE D-1. 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 D-2. ILLUSTRATION OF THE PROBABILITY THAT f(t) OCCURS
IN THE INTERVAL fn- TO f-j + 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
(D-2)
where the summation is over all intervals f.- in the range -°° < f < °°.
J
Let the mean value of the random sequence f(t-j), i = 1,2,... ,m, generated
above be noted by
(D-3)
Then by definition the ensemble mean value of the random variable f(t-j) is
f = lim f = lim 2,
" Af+0 j
(D-4)
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. (D-l) and (317),
we obtain
f = lim 2, f-;
Af+0 j J
m(e.
T
Comparing this with Eq. (D-2), we see finally that
f
,+T
f(t) dt = T li
m
n
1 2, f(ti
(D-5)
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(t-j) formed by picking points t-j
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 one-third. 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 one-dimensional tubular reactor having a
head made of some 100 small nozzles (Figure E-l). Reactants are fed through
alternate jets to simulate a cross-sectionally 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 E-l. 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:
(E-2)
2
(E-3)
Since the species are not premixed,
= - . (E-4)
Combining Eqs. (E-l) through (E-4), we obtain
-------
221
(A B (X)> (E-5)
ii' — * x / v x »•«••* \ / w \ ^ /
-------
222
o
•
-------
223
VERY RRPID RERCTION
STGICHIOMETRIC FEED
CD MRO, RE=1460
— CURVE FITTING
0.0
FIGURE E-3. DECAY OF CONCENTRATION WITH DISTANCE
IN MAO'S REACTOR
-------
224
VERY RRPIO RERCTION
STOICHIOMETRIC FEED
O MRO, RE=2i430
— CURVE FITTING
0.0
FIGURE E-4. DECAY OF CONCENTRATION WITH DISTANCE
IN MAO'S REACTOR
-------
225
APPENDIX F
ESTIMATION OF THE ORDER OF MAGNITUDE OF
THE TIME SCALE OF DECAY OF A"2
-------
226
APPENDIX F
ESTIMATION OF THE ORDER OF MAGNITUDE OF
THE TIME SCALE OF DECAY OF A"2
When the first term on the righthand side of Eq. (156) is expressed
in the form (157), we find that, in the absence of any mean transport,
chemical decay, or sources of subgrid-scale concentration fluctuations,
iV*"*)
an initially present mean square fluctuation field A" decays according to
A" , (F-l)
where A is some length scale and KA is the eddy diffusivity. In other
words, the time scale of decay of A"? due to turbulent diffusion is of
the order of
(F-2)
We wish to show from an estimate of T that A is of the order of the size
-^
of the grid mesh upon which the mean concentration field A is defined.
For simplicity we consider a one-dimensional problem governed by
the equation
= K^- + «(z)6(t) (F-3)
o Z
with initial and boundary conditions
c(z,0) = 0 (F-4a)
lim c(z,t) = 0 , (F-4b)
-------
227
where K is the turbulent diffusivity, assumed to be a constant, and 5 is
the delta function. The solution of (F-3) and (F-4a,b) is found to be
c(z,t) = (4TiKt)"1/2 expf- fiTT- ) .
V KV
Now let
Z+AZ
c(z,t) =
2Az /
c(z',t)dz'
Z-AZ
Averaging (F-5) in this manner we obtain
(F-5)
(F-6)
c(z,t) =
4Az
erf ^^ -
z-Az
(F-7)
where
erf(x) =
/
-^
e dx
is the standard error function. Note that Eq. (F-7) is the solution
of the equation obtained by averaging Eq. (F-3), that is
(F-8)
where
= 1, -Az <_ z < Az;
;
0, otherwise.
-------
228
Let us assume that Eq. (F-8) is a grid model representation of Eq. (F-3)
and that c, as given by Eq. (F-7), is the model's output. The subgrid-
scale variations are therefore [from Eqs. (F-5) and (F-7)]
c"(z,t) = -==Lr exp j- — .
. V4i:Kt V 4Kt/ 4Az
erf f£y£
_ era
er9 .——
• (F-9)
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. Moreover,
for fixed Az the SSV in this example decrease with time.
To obtain a more concise description of the relative magnitude of the
subgrid-scale 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
- c" o.t
= „.,
c 0,t
(F-1Q)
From the expressions for c" and c we obtain
where
a*erf f -TT
- 1
(F-ll)
T* =
(F-12)
Az
In words, a* is the approximate half-width of the pollutant cloud, normalized
by the grid mesh size AZ, at time t; It is evident from Eq. (F-ll) that a*
-------
229
is the key parameter which determines the relative size of the SSV in this
instance. On evaluating Eq. (F-ll) for several values of a*, we obtain
p = 10 when 0* = 0.1; p = 0.34 when 0* - 1; and P = 0.08 when a* = 2. Thus,
the time scale of decay of p is of the order
T
P 4K K ^r~'°
Now since c (0,t), which enters in the definition of p [in Eq. (F-10)J, also
decays with time, the decay rate of c"(0,t) can be no faster than that of p.
Furthermore, since c"(0,t) represents the maximum amplitude of c^(z,t), it
2
is not unreasonable to assume that the time scale of decay of c" is com-
parable to that of c"(0,t). From all of these considerations we conclude
that
and hence that the length scale A, defined in the beginning of this appendix,
is of the order of AZ.
-------
230
REFERENCES
Batchelor, G. K. (1967), "Small Scale Variation of Convected Quantities Like
Temperature in a Turbulent Fuel, Part I," J. of Fluid t-lech. , Vol. 5, p. 113.
(1950), "Application of the Similarity Theory of Turbulence to Atmos-
pheric Diffusion," Quart. J. of the Royal Meteor. Soc., Vol. 76, pp. 133-146.
Bendat, J. S., and A. G. Piersol (1968), "Measurement and Analysis of Random Data,"
John Wiley and Sons, New York, New York.
Briggs, G. A. (1968), in Meteorology and Atomic Energy, David Slade, ed.
Chen, W. H., and J. H. Seinfeld (1974), "Estimation of Spatially Varying Parameters
in Partial Differential Equations," Int. J. Control, Vol. 15, pp. 487-495.
Deardorff, J. (1972), "Numerical Investigation of Neutral and Unstable Planetary
Boundary Layers," J. Atmos. Sci., Vol. 29, pp. 91-115.
(1970), "A Three-Dimensional Numerical Investigation of the Idealized
Planetary Boundary Layer," Geophys. F1uid Dyn. , Vol. 1, pp. 377-410.
Donaldson, C. duP., et al. (1973), "Atmospheric Turbulence and the Dispersal of
Atmospheric Pollutants," Report EPA-R4-73-016a-c, Environmental Protection
Agency, Research Triangle Park, North Carolina.
Donaldson, C. duP. (1969), J. of AIAA, Vol. 7, pp. 271-278.
Golder, D. (1972), "Relations Among Stability Parameters in the Surface Layer,"
Boundary Layer Meteorology, Vol. 3, pp. 47-58.
Hanjalic, K., and B. E. Launder (1972), "A Reynolds Stress Model of Turbulence
and Its Application to Thin Shear Flows," J. Fluid Mech. , Vol. 52, pp. 609-638
Hilst, G. R., et al. (1973), "A Coupled Two-Dimensional Diffusion and Chemistry
Model for Turbulent and Inhomogeneously Mixed Reaction Systems," Report
RA-73-016C, Environmental Protection Agency, Research Triangle Park, North
Carol ina.
Hinze, J. 0. (1959), Turbulence, McGraw-Hill Book Company, New York, New York.
Kreith, F. (1966), Principles of Heat Transfer, Second Edition (International
Textbook Company, Scranton, Pennsylvania).
Lamb, R. G. (1974), "Lagrangian Analysis of a Passive Scalar in Turbulent Fluid,
Part 1: General Theory," J. Fluid Mech., submitted for publication.
(1971), "Numerical Modeling of Urban Air Pollution" PhD. dissertation.
Department of Meteorology, University of California, Los Angeles, California.
-------
Lamb, R. G. (1973), "Note on the Application of K-Theory to Diffusion Problems
Involving Nonlinear Chemical Reactions," Atmos. Env. , Vol. 7, pp. 257-263.
Lamb, R. G., and J. H. Seinfeld (1973), "Mathematical Modeling of Urban Air
Pollution — General Theory," Environ. Sci. Techno]^, Vol. 7, pp. 253-261.
Lamb, R. G., and M. Neiburger (1971), "An Interim Version of a Generalized
Urban Air Pollution Model," Atmos. Environ. , Vol. 5, pp. 239-264.
Lettau, H. H. (1959), "Wind Profile, Surface Stress, and Geostrophic Drag Coef-
ficients in the Atmospheric Surface Layer," Advances in Geophysics, Vol. 6,
pp. 241-256.
Lions, J. L. (1971), Optimal Contro1_of Systems Governed by Partial Differential
Equations (Sprinqer-Verlaq, New York, New York) 396 pages.
Lumley, J. L. (1970), "Toward a Turbulent Constitutive Relation," J. Fluid Mech.,
Vol. 41, pp. 413-434.
Lumley, J. L., and B. Khajeh-Nouri (1973), "Computational Modeling of Turbulent
Transport," Second IUTAM-IUGG Symposium on Turbulent Diffusion and Environ-
mental Pollution, Charlottesville, Virginia.
Mao, K. (1969),' "Chemical Reactions with Turbulent Mixing," Ph.D. thesis,
Carnegie-Mellon University, Pittsburg, Pennsylvania (173 pages).
McElroy, J. L. (1969), "A Comparative Study of Urban and Rural Dispersion,"
J. Appl. Meteor., Vol. 8, pp. 19-31.
McKelvey, K. N. (1968), "Turbulent Mixing with Chemical Reactions," Ph.D. thesis,
Ohio State University, Columbus, Ohio (123 pages).
Monin, A.^S., and A. M. Yaglom (1971), Statistical Fluid Mechanics, MIT Press,
Cambridge, Massachusetts.
O'Brien, E. E. (1968), "Closure for Stochastically Distributed Second-Order
Reactants," Physics of Fluids, Vol. 11, p. 1883.
Orszag, S., and M. Israeli (1974), "Numerical Simulation of Viscous Incompressible
Flows," Ann. Rev, of Fluid Mech., Vol. 6, p. 281.
Ragland, K. W. (1973), "Multiple Box Model for Dispersion of Air Pollutants from
Area Sources," Atmos. Environ. , Vol. 7, pp. 1017-1032.
Shir, C. C., and L. J. Shieh (1973), A Generalized Urban Air Pollution Model and
Its Application to the Study of S02 Distributions in the St. Louis Metropolitan
Area," Report RJ 1227, IBM, San Jose, California.
Shu, W. R. (1975), Ph.D. dissertation, in progress.
-------
232
Toor, H. L. (1969) Ind. Eng. Chem. Fundam., Vol. 8, pp. 655-659.
(1962), "Mass Transfer in Dilute Turbulent and Nonturbulent
Systems with Rapid Irreversible Reactions and Equal Diffusivities,"
A^LCh.E. Journal, Vol. 8', pp. 70-78.
Turner, D. B. (1969), "Workbook of Atmospheric Dispersion Estimates," Publi-
cation AP-26, Environmental Protection Agency, Research Triangle Park,
North Carolina.
Vassilatos, G., and H. L. Toor (1965), "Second-Order Chemical Reaction in a
Nonhomogeneous Turbulent Fluid," A.I.Ch.E. Journal, Vol. 11, pp. 666-673.
Zakanycz, S. (1971), "Turbulence and the Mixing of Binary Gases," Ph.D. thesis,
Ohio State University, Columbus, Ohio (172 pages).
-------
233
TECHNICAL REPORT DATA
(f'lcasc read Inunii'liun^ on l/ic /i i im1 be tore coinpli InifJ
.'.-'iOO/.'4-76rQ.16_c_
. , ,-.'. ^ bun 1 II LE
CONTINUED RESEARCH IN MESOSCALE AIR
il.I.KIlON SIMULATION MODELING. VOLUME III. Modeling
!' v'i i-roscnle Phenomena
5. REPORT DATE
_ May T9 7 6
6. PERFORMING ORGANIZATION CODE
R. G. Lamb
8. PERFORMING ORGANIZATION REPORT NO.
EF75-25
>•(. ill OH VI.NT, ORG \NIZATION NAME AND ADDRESS
SYSTEMS APPLICATIONS, INC.
950 NORTHGATE DRIVE
SAN RAFAEL, CALIFORNIA 94903
I/ ;,»0
-------