-------
39
y,
O.I
0.01
i i .1 i
O
A
o
A
=0.55
ugO.5,' ft= 1000:
2000
= 0.35
?0.44,h=u*/f
- McElroy- Pooler cry h = 1000 meter
neutral and unstable
o Deardorff ' o-x fi/L=-4.5 (longitudinal o-x)
Deardorff crx
i i i i
= 0 (longitudinal crx )
i i i i
0.01
O.I
t
FIGURE 12. COMPARISON OF THE COMPUTED NONDIMENSIONALIZED STREAMWISE
ROOT-MEAN-SQUARE PARTICLE DISPLACEMENT ax (CIRCLES AND
TRIANGLES) WITH THE MCELROY-POOLER LATERAL DISPERSION
PARAMETER ay FOR NEUTRAL AND UNSTABLE CONDITIONS AS A
FUNCTION OF TRAVEL TIME t
-------
40
that since the empirical a's are in close agreement with the numerically
calculated values of x2'5, z2'1, and so forth, major errors in the puff
model predictions reflect departures of the turbulence characteristics
from those consistent with the assumed Gaussian form of p. Similar ar-
guments apply to the plume formula. Let us examine these errors.
For this purpose, we set up the standard plume model using Pasquill-
Gifford dispersion data, and the puff model of Lamb and Neiburger (1971)
using the McEl roy-Pooler data, to simulate the same two-dimensional problem
that was treated earlier in Section C. In each case, we nondimensionalized
the empirical data in the manner portrayed in Figures 10, 11, and 12 so
that each set was in closest agreement with the computed particle statis-
~Z
tics x , z , and so forth. We should point out that the plume model does
not account rigorously for the effects of the inversion layer present in
this problem. However, in the trial calculations presented here, this de-
ficiency is not serious, because within the range of downwind distances
treated, the inversion layer has little overall influence on the concen-
tration distributions.
Figure 13 presents the results of the calculations for the neutral
case, and Figure 14 gives those for the unstable case. As in the previous
figures, we have plotted the model predictions in terms of their fractional
departure e, Eq. (42), from the corresponding numeri co-empirical solutions
of Eq. (16). The latter are given in Figures 8(a) and 9(a). We consider
only the puff model results here, because these results provide a better
assessment of the accuracy of the assumed Gaussian form of p due to the
closer agreement between the McEl roy-Pooler data and the computed particle
statistics.
As shown in Figure 11, for the neutral case the empirical a and the
—^— i L-
calculated z2"5 profiles are nearly superposed for travel times t < 0.1,
which corresponds to distances x from the source in the range 0 < x < 4.
-------
II
u*
(a) Nondimensionalized Cross-Wind Integrated Concentrations
Predicted by the Puff Model, Expressed as the Percentage
Deviation (i.e., lOOe) from Those Shown in Figure 8(a)
Which Were Obtained from Eq. (16) for the Case of Neutral
Stability
0.4- -
(b) Results of the Plume Model Calculations Expressed as in
Part (a)
FIGURE 13. COMPARISON OF THE PREDICTIONS MADE BY THE PUFF MODEL AND
THE PLUME MODEL FOR THE CASE OF NEUTRAL STABILITY. The"
puff model used neutral, travel-time-dependent McElroy-
Pooler data; the plume model used Pasquill-Gifford Class C
data and an expanded x-axis (see the text).
-------
0
16
(a) Nondimensionalized Cross-Hind Integrated Concentrations
Predicted by the Puff Model, Expressed as the Percentage
Deviation (i.e., 100 ) from Those Shown in Figure 9(a),
Which Were Obtained from Eq. (16) for the Unstable Case
(b) Results of the Plume Model Calculations Expressed as in
Part (a)
FIGURE 14. COMPARISON OF THE PREDICTIONS MADE BY THE PUFF MODEL AND
THE PLUME MODEL FOR THE UNSTABLE CASE. , The puff model
used "slightly unstable" McElroy-Pooler data; the plume
model used Pasquill-Gifford Class B data and an expanded
x-axis (see the text).
-------
43
However, Figure 13(a) reveals that in this range the predictions of the
puff model are considerably in error. Concentrations are underestimated
on the centerline of the plume and are greatly overestimated at the
edges, especially along the plume edge that touches the ground. This
behavior suggests that the diffusing particles behave as if a weak re-
storing force were preventing them from wandering too far away from their
level of release, i.e., the plume axis. Thus, particles that have wan-
dered a distance z^2, say, from the plume axis have a larger probability
of moving back toward the plume centerline than of continuing on to still
more remote points. In contrast, the Gaussian density implies that at
each instant the particles have equal probabilities of moving toward or
away from the plume axis. In short, the Gaussian density has a much
larger kurtosis or flatness factor than the true probability density p
does under neutral conditions.
Turning next to the unstable case, we see in Figure 11 that the measured
a values and those computed for z" are in excellent agreement at all travel
times (and distances from the source), except those at the extreme upper limit
of the range treated. Figure 14(a) reveals that the errors in the puff model
predictions are smaller in this instance than those in the neutral case but
that there is still a tendency for the predicted concentrations along the
plume centerline to be too small and for those at the edges to be too large.
One factor contributing to these errors is the assumption in the puff model
that the joint moment xz~ is zero. Actually, this moment is strongly positive;
because of the rapid increase in wind speed above the level of release (namely,
s*
Zpp. = 0.025h) of the particles in the unstable case (see Figure 5), particles
that are displaced upward are also systematically displaced downstream. Con-
sequently, at any distance x from the source, a significant fraction of the
particles that normally would be found at the plume edge are found farther
downstream, and in their places are particles from points upstream of x, which
thus have smaller vertical displacements. As a result, the effective width of
~2)^
the plume i.s consistently smaller than one would predict on the basis of z 2
alone. This phenomenon is not as pronounced in the neutral case, because the
wind shear at the level of release (znri = O.lu*/f) of the particles is much
smaller (see Figure 5).
-------
44
The analyses just completed illustrate the inadequacy of the basic
assumption underlying the plume and puff models that the probability
density p entering in Eq. (16) is of a Gaussian form. However, despite
these limitations, the plume and puff models remain attractive formula-
tions in applied studies because of their mathematical simplicity. For
this reason, it is useful to inquire whether profiles exist that will
bring the errors in the predictions of these two models within the range
of acceptability.
To explore this question, we used the variances a (T) and a (T) as
A. L-
adjustable parameters for obtaining a least-squares best fit of the
Gaussian density function to the numerico-empirical expression p. Figures
15 and 16 present the resulting "optimal" profiles of a and a , respec-
X £-
tively. In each case, the corresponding empirical data of McElroy and
Pooler, presented earlier in Figures 11 and 12, are shown for comparison.
Since the method used to obtain the empirical data is similar to that used
here to.obtain the optimal a and a values, it is surprising that the lat-
A L.
ter do not compare as well with the empirical data as do the computed root-
mean-square particle displacements shown in Figures 11 and 12.
In a similar manner, we adjusted the dispersion parameter a (x) to ob-
tain the least-squares best fit of the plume formula to the numerico-empirical
concentration distributions given in Figures 8(a) and 9(a). Figure 17 pre-
sents the results, along with the empirical Pasquill-Gifford data for compar-
ison. We subsequently inserted these "optimal" dispersion coefficients into
the puff and plume models and performed new sets of calculations to compare
with the numerico-empirical concentration profiles given in Figures 8(a) and
9(a). The results, given in terms of the fractional error e used earlier,
are presented in Figure 18 for the neutral case and in Figure 19 for the un-
stable case. Comparing Figures 13(a) and 14(a) with Figures 18(a) and 19(a),
we find that although the optimal a and a profiles improved, the overall
A i-
performance of the puff model, errors in the predicted values remain intol-
erably large in certain regions. This is especially true for the ground-level
concentrations predicted in the neutral case [see Figure 18(a)].
Figures 13(b) and 18(b) reveal that using the "optimal" az(x) profiles,
the plume model fails by a considerable margin to achieve the level of accur-
acy that was attained in the earlier calculations. The reader will recall
that in the previous computations, we expanded the x-axis to account for
-------
o Optimal h/L = -4.5
Optimal h/L = 0
I8
i= ±O.OlJ
0.05-
FIGURE 15. OPTIMAL az(t) PROFILES FOR USE IN THE PUFF MODEL
COMPARED WITH THE EMPIRICAL DATA OF MCELROY AND POOLER
-------
46
0.5
O.I
u*=0.5,h=IOOO
u*=0.44th=u*/f
o Optimal h/L=-4.5
A Optimal h/L=0
McElroy-Pooler
(Neutral S Unstable)
0.5
0.
0.5
FIGURE 16. OPTIMAL STREAMWISE DISPERSION ax(t) COMPARED WITH
THE EMPIRICAL LATERAL DISPERSION DATA OF MCELROY AND POOLER
-------
47
0.5
cr2 O.I
0.05
h=IOOO
= 0.5
o Optima! h/L - -4.5
A Optimal h/L=0
---- GLASS B
- CLASS
-- CLASS
oj
Gifford
_J I
10
FIGURE 17. OPTIMAL az(x) PROFILES FOR USE IN THE PLUME
FORMULA COMPARED WITH THE EMPIRICAL DATA
OF PASQUILL AND GIFFORD
-------
4'B
0.4
0.3
0.2
O.I
f- -0.10
-51/?
e2 = 0.77
8 10 12 14
16
(a) Results of the Puff Model
i=0.02
il
u*
FIGURE 18.
(b) Results of the Plume Model
THE FRACTIONAL ERROR (lOOe) RESULTING IN THE PUFF MODEL AND
THE PLUME MODEL USING THE OPTIMAL a PROFILES SHOWN IN FIGURES
15 THROUGH 17 FOR THE NEUTRAL CASE. Errors are relative to
numeri co-empirical concentration distribution given in
Figure 8(a).
-------
49
A
~
h
0.8-
0.2-
8 10 12
X
(a) Results of the Puff Model
14
(b) Results of the Plume Model
16
FIGURE 19,
THE FRACTIONAL ERROR (lOOe) RESULTING IN THE PUFF MODEL AND
THE PLUME MODEL USING THE OPTIMAL a PROFILES SHOWN IN FIGURES
15 THROUGH 17 FOR THE UNSTABLE CASE. Errors are relative to
the numerico-empirical concentration distribution given in
Figure 9(a).
-------
apparent differences in the mean particle velocities implicit in the empirical
data and in the numerical turbulence model. We also made arbitrary adjust-
ments of the atmospheric stability. In contrast, we did not employ any of
these artifices in computing the optimal a (x) values. These facts, together
with the error values given in Figure 18(b) suggest that the basic plume
formula does not provide a wholly adequate description of atmospheric
diffusion, at least under neutral stability conditions. However, under
unstable conditions, the accuracy is better, as indicated by Figure 19(b).
E. IMPLEMENTATION OF THE OPTIMAL DIFFUSIVITY AND
WIND SHEAR PROFILES IN AN AIRSHED MODEL
Since the optimal diffusivity profiles derived in Section C are based
on a rather small ensemble of particle trajectories, the statistical sig-
nificance of these profiles is limited. In addition, the question arises
of whether the profiles vary significantly with changes in the source
height. (In our studies, data from only one release height were available
for each stability case.) The data needed to resolve these problems are
currently being collected using Deardorff's newest model, and we plan to
continue studies similar to those described above but with a much wider
scope.
For the present, however, we are developing the methods necessary to
implement these diffusivity and wind shear profiles in operational airshed
models. For this purpose, we have fitted fourth-order polynomials to the
wind shear and optimal diffusivity profiles derived earlier. Table 2 pre-
sents these in nondimensional form. The only task necessary to use these
results in a diffusion model is to convert the profiles into dimensional
forms. Thus, denoting dimensional variables by the caret ("), we have [cf
Eq. (22)]
K (--ju z, , unstable case (z./L = -4.5) ; _ (43)
Z \ i- */ f I 1
Kz(z) = <
neutral case (z./L = 0)
-------
51
Here KZ(Z) is the nondimensional diffusi vity, given in Table 2, at height
z (nondimensional) and u , z., f, and L are defined as before. Similarly,
"ft \
we have
(44)
where u is the friction velocity, u is the wind speed profile given
"K ^
in Table 2, and h = z. in the unstable case and u /f under neutral condi-
tions. Note that the wind profile u given in the table is for a particu-
A ^ A.
lar value of R = h/z~, where z~ is the surface roughness. For different
values of R, the profiles must be adjusted as described below.
In air pollution modeling studies, one normally has available only
*
z., f, z~, and u,n, where the last term denotes the wind speed at 10 m .
Atmospheric stability data are seldom available from measurements. Rather,
Pasquill's technique of categorizing the stability into Classes A through F
using ground level wind speed and insolation observations is frequently em —
ployed. Recently, Golder (1972) obtained approximate relationships between
Pasquill 's stability classes and the Monin-Obukhov length L. Thus, we
assume that estimates of z./L can be obtained using Golder's formulas and
Pasquill's stability classes. The friction velocity u^ must be determined
before the wind and diffusivity profiles can be converted into dimensional
forms .
An often used formula for u is
= CDug
where Cn is a drag coefficient for the boundary layer and u is the geostrophic
wind speed. Lettau (1959) proposed the formula
Ii u \ ]-l
Io9in f ~1-8 ' neutral case (46)
IUVZ0/ J
A technique for estimating the surface roughness z~ has been proposed by
Lettau (1970).
-------
Table 2
OPTIMAL PARAMETERS FOR USE IN THE DIFFUSION EQUATION AND IN THE
GAUSSIAN PUFF AND PLUME MODELS (NONDIMENSIONAL VARIABLESt)
Model
Gaussian
puff
Parameter
Stability
Neutral - = 0
1.5 x 10
0.0045 + 0.4r, 0.022 < t <_ 0.079
0.026 + 0.129-t, 0.079 < T < 0.32
Unstable j1 = -4.5 —- = 6.8 x 106
. I; zo •
= -0.0023 + 0.65t - 1.89r2 + 5.32T3,
0.04 < T < 0.52
«X(T)
1.85T, 0 < T <_ 0.045
0.045 + 0.92-r, 0.045 < T < 0.32
- 0.045 + 1.307-r, 0.05 < T < 0.7
Gaussian
plume
oz(x)
0.032 + 0.0044X, 0.23 < x < 4.0
= 0.075 + 0.0167x, 0.5 < x < 5.5
- -0.065 + 0.0435X, 5.5 < x < 9.5
Diffusion
equation
« 7.396 x 10"4 + 6.082 x 10"Zz
+ 2.532z2 - 12.72z3 + 15.T7z4,
0 < z < 0.45;
-6.934-x 10"3 + 0.6113Z + 3.297z2
-6.442z3 + 3.153z4,
0.02 < z < 1.0
0, z > 0.45.
u(z)-
29.82 + 213.2z - 1989z'
26.22 + 153.2z - 1428z':
8743z3 - 14670z4, 0.022 < z < 0.21; + 5541z3 - 7523_z4,_0.025 < z < O.J3
u(0.21), z > 0.21.
u(0.3), z > 0.3.
t All lengths were made nondimensional by zi » Inversion height in the unstable case and by u^f in the
neutral case. Times were made nondimensional by Zj/u* (unstable) and 1/f (neutral).
5 The polynomial given for the neutral case is only qualitatively accurate. See Appendix A for a more
precise formulation.
en
-------
53
for neutral conditions. Our work with Deardorffs boundary layer model has
led to a confirmation of this expression. We have also found that in the
unstable case, the drag coefficient is given approximately by
CD = 0.156
/ i \ I"1
loglo(z )~ 1>18 , unstable case . (47)
Because the geostrophic wind speed u is often not known in air pollution
studies, Eq. (45) is difficult to use in diffusion modeling applications.
Consequently, we attempt below to obtain u from the 10-m wind u,n and the
wind speed profiles that we have available.
Consider first the unstable case z./L = -4.5, which is representative
of the mixed layer over urban regions during much of the daylight period.
For the unstable case, we have from Table 2,
~/L\
u zi ?
~^-LL = u(z) = 26.22 + 153.2z - 1428z^
*
+ 5541z3 - 7523z4
(48)
0.025 < z/z. < 0.3
= 0(0.3) , z = z/z. > 0.3
Note that this profile is valid only for z^ 0.025, R = Z./ZQ = 6.8 x 106, and
z./L = -4.5. For different values of R, the wind speed profile U(Z,ZQ) is (for
a given ZQ)
u(z,zn) = u(z) - p£n(6.8 x 106 zn) , 0.025 < z < 0.3 , (49)
U K V U/ ZQ < .004 ;
where u(z) is given by Eq. (48) and k is the von Karman constant (-,35). (Note
the restriction on the size of Z.)
Let
Z10 = (50)
i
-------
represent the nondimensionalized 10-m anemometer height,
from Eq. (49)
Then we have
u(z,n,zn) =
J10
10' 0' u
Therefore,
u.
u = .,10
* D(zio'zo)
(51)
Note that the expression in Eq. (49) for u(z,zn) is valid only for heights
z > 0.025. If z,~ is smaller than this value, we are in the surface layer,
where similarity theory gives (see Table 1)
u(z,zn) = |- 2/tan"1 x - tan^xj + £n f;
u K|_ \ u/ yi
ZQ < z < 0.025
where
y —
,1/4
- 15(z
(52)
1/4
(53)
In summary, we obtain the friction velocity u^ in unstable cases from the
formula
"10
"* = "^^O^O7
where u,n is the wind speed at the anemometer height and u is given by
Eqs. (48) and (49) for z]Q in the range 0.025 < ZIQ < 0.3 and by Eqs. (52)
and (53) when zn < z,Q < 0.025. Note that the surface roughness z^, the
mixing height z., and the Monin-Obukhov length L must be known to obtain u
-------
55
In the neutral case (z./L = 0), the evaluation of u^ is somewhat
simpler because the anemometer height is always within the surface layer.
Here similarity theory gives (see Table 1)
1 /z + zn\ /z.
u(z,zn) = r £n - - , h-L = OJ neutral stability
\ zr> / \L /
\ 0 / \ / (54)
Thus, if Z-|Q = 10 meters and u,Q is the 10-m wind speed, we can assume
that
> neutral stability . (55)
It is important to note that Eq. (54) is very accurate up to a height
z * 0.03 u^/f. Since f « 10~4 sec"1 for latitudes within the United
States and u > 0.1 m/sec over urban regions during most of the daylight
hours, it is easy to see that 10 m is a value that is well within the
range of validity of Eq. (54).
Having established formulas for the friction velocity u^, we can con-
vert the nondimensional wind speed and diffusivity profiles of Table 2 into
physical forms quite easily using the computer. Our approach is to define
FORTRAN functions USTAR, DKZ, and UBAR, which calculate the friction velocity,
diffusivity, and wind speed, respectively, given the following data:
uin = anemometer wind speed at height z,,.,
zn = local surface roughness,
z. = local inversion height,
z./L = stability parameter (recall that the Monin-Obukhov
1 length L is assumed to be known),
z = level at which the diffusivity or wind speed is required,
f = 2fi sin = Coriolis parameter.
-------
56
In the execution of the pollution model, the function USTAR is called
for each point on the horizontal plane to compute the local friction veloc-
ity. Thus, we have
u^ = USTAR(U10,ZO,ZIOVL,ZI,Z10)
where the arguments of USTAR are given the proper local values. Once u is
determined, the functions DKZ and UBAR can be referenced as follows:
K (z) = DKZ(Z,ZI,USTAR,F,ZIOVL)
u(z) = UBAR (Z,ZO,ZI,ZIOVL,F,USTAR)
These functions merely use their arguments to translate the polynomials
given in Table 2 into dimensional values. Note that the functions DKZ and
UBAR simply replace the diffusivity and wind speed variables in the program.
Listings of the FORTRAN programs of each of these three functions are
provided in Appendix A. Note in the listing of UBAR that in the portion of
the surface layer not covered by the profiles given in Table 2, we use well-
known analytical expressions for u(z). Equation (52) is an example of such
an expression. Further information in this regard can be obtained from
Ragland (1973).
At present, the three functions USTAR, DKZ, and UBAR are applicable only
to the stability cases z./L = 0 (neutral) and z./L = -4.5 (slightly unstable).
The data on which the profiles of Table 2 are based are currently restricted
to these conditions. We hope to be able to extend the scope and accuracy of
this work in the future. Deardorff (1972) has found, for example, that in all
situations where -z./L is greater than about 5, the proper velocity scale is
-------
57
rather than u . If this is true, then one K profile should suffice for
"rT i-
all unstable conditions in the range z./L < -5, but this remains to be
checked. Also, our K profiles were obtained for fixed source heights.
In principle, the diffusivity should be independent of source height;
but because of the heuristic nature of the concept of a turbulent diffu-
sivity, this may not be the case. Further work is therefore needed to
ascertain the influence of source height on the optimal K profiles.
-------
58
III DEVELOPMENT OF A CLOSURE APPROXIMATION
FOR THE CONCENTRATION FLUCTUATION TERMS
IN THE GOVERNING EQUATIONS
A. INTRODUCTION
Snapshots of smoke plumes in the atmosphere always reveal an erratic
distribution of smoke concentration within the local vicinity of the smoke
source. In contrast, time exposures usually reveal a smooth distribution
that closely resembles the classical Gaussian profile. If two such photo-
graphs of the same plume were superposed, the difference in apparent smoke
concentration at any point would represent the instantaneous magnitude at
that point of the so-called turbulent concentration fluctuation. That is,
if A denotes the instantaneous concentration, portrayed by the snapshot, and
A denotes the mean, revealed by the time exposure, then
A1 = A - A (56)
is the concentration fluctuation.
In situations where nonlinear chemical reactions occur among consti-
tuents of the plume or between the plume's constituents and chemical species
entrained from the ambient atmosphere, the presence of concentration fluctua-
tions can have a profound impact on the behavior of the mean concentrations
of all the reacting species. This can be demonstrated mathematically by the
following simple example. Suppose that two species, A and B, undergo the
bimolecular reaction
A + B
-------
59
with rate constant k. The equation governing the instantaneous concentration
of A is, therefore,
{£= -kAB . (57)
If the two species are not homogeneously distributed, local fluctuations
A1 and B1 exist. The mean concentration A" is then governed by the
equation
__ _
= -kAB - kA'B1 . (58)
This equation follows from Eqs. (57) and (56) and the assumption that
A1 = B1 = 0. Equation (58) indicates that, depending on their amplitude,
the concentration fluctuations can actually dominate the evolution of the
mean concentration distribution.
Since the term A'B' is an additional unknown variable, Eq. (58) cannot
be solved until another expression relating A'B' to A~ and B" or to some other
known quantities is found. This is known as the closure problem. The
search for such relationships, or closure approximations as they are generally
called, has been undertaken by several investigators. Of these, the most
sophisticated work has been done by O'Brien (1968) and Hilst et al.
(1973). Both of these studies, which are similar in nature, follow along
the lines of the so-called invariant modeling approach. We outline this
method below, using the work of Hilst et al. as an example.
In the context of the problem represented by Eq. (58), the first step
is to formulate the equation governing A'B' and to add it to the system of
equations governing A and B. We then have [repeating Eq. (58)]
^ = -kAB - kAW (59a).
^:= -kAB - kFB^ (59b)
-------
'—- = -k|_AB'2 + BA'B1 + BA'2 + AA'B'
(59c)
+ A'B'2 + B'A'2J
2 2
The last equation involves the mean squares A' and B1 , and so their
governing equations must also be considered:
'2
dA'
dt
= -2k[BA'2 + AA'B' + A'2B'J . (59d)
21- . 1
„„ . o" 9"
^-= -2k[AB^ + BA'B1 + A'B'^J . (59e)
2
The presence of third-order moments, such as A'B1 , in Eqs. (59c) through
(59e) precludes a solution of the system of equations represented by
Eqs. (59a) through (59e). Note that the problem here is identical to the
one we faced with Eq. (58), i.e., too many unknowns for the available set
of equations. It is apparent, therefore, that merely writing down the
governing equation for each new unknown that appears does not in itself
resolve the difficulty, for this process produces an infinite hierarchy
of equations. At some point, a hypothetical expression must be introduced
to circumvent the hierarchy and to produce a closed, finite set of equations.
It is at the point represented by the system [Eqs. (59a) through (59e)] that
Hilst et al. introduced their closure hypothesis in the form
2 '2 M
l2R, _ AB , , A
-
-------
where
M =
0, if (A'2/A2) (B'2/B2) < 1
1, otherwise
Equations (59f) and (59g) are not based on any physical concepts or
on experimental observations. Rather, they are simply mathematical approxi-
mations that are designed to satisfy the following realizability conditions,
= 0
if A'B1 = -AB
A2B = A2B , if A'2 = 0 ;
0
2B = A2B + A|2B +
AB = AB + A|B + 2A'B'A
if A'2B' = 0 ;
(60)
We should note that Eqs. (59f) and (59g) are not the only relationships
that satisfy Conditions (60), nor are these conditions the only ones that
the statistics of A and B satisfy.
Although Hilst et al. have found good agreement between the pre-
dictions of their model [Eq. (59)] and the corresponding analytic solu-
tions of certain simplified chemical reactors, from the standpoint of
urban pollution modeling, their approach is not altogether acceptable. The
objection is that in a system involving m nonlinear reactions among n species,
their technique introduced n + m additional differential equations into the
existing set of n to be solved. Thus, in a typical case where there are
five pollutants and five reactions (in which concentration fluctuations
are important), the number of differential equations to be solved is 15
instead of 5. Needless to say, an increase in computational effort of
this magnitude cannot be accommodated in urban diffusion models given pre-
sent computer speed.
-------
62
Since none of the existing closure approximations for the turbu-
lent concentration fluctuation terms possess both the degree of accuracy
and the ease of implementation that are required in pollution simulation
studies, we set for ourselves the task of developing a new approximation
that would more closely satisfy these needs. Our primary constraint was
that the new scheme should involve no additional differential equations.
The existence of such a scheme has been suggested by our work in the pre-
vious chapter with the analogous term u]c'. In Chapter II, we found
that by a proper choice of the functional form of the diffusivity KZ, the
approximation
could be made sufficiently accurate to be useful in diffusion modeling
applications. Moreover, this formulation does not require that additional
equations be solved. Using this work as a guide, we set out to find an
expression of the form
A'B1 = f(A,B,known parameters)
We describe the expression that we have developed and our attempts to
demonstrate its accuracy in the remainder of this chapter.
B. DERIVATION OF AN APPROXIMATION FOR CONCENTRATION FLUCTUATION TERMS
SUCH AS A'B1
The following simple example illustrates both the reasoning behind
our approach and the essence of the closure problem itself.
Suppose that rectangular pulses or slugs of the species A and B are
present in a one-dimensional system as shown in Figure 20. The mean con-
centration A of A is defined by
/
XO+L
xo
-------
£3
CONCENTRATION
FIGURE 20. TWO INITIAL SLUGS OF SPECIES A AND B
IN A ONE-DIMENSIONAL VESSEL
and a similar definition applies for B. These results follow immediately
from Figure 20. (Although in this example L represents an averaging
length interval, it could just as well denote time or ensemble averaging
domains because the results that we derive below are the same in all cases.)
Since A and B are not uniformly distributed within L, concentration fluctua-
tions exist. At points outside the slug of A, the fluctuation A' has a
magnitude
A1 = A - A = -
VA
(62)
and at points within the slug,
A' •
(63)
Similar expressions pertain to the fluctuations B1
Suppose that both A and B undergo first-order decay reactions whose
rates are given by
dA
dt
dB
dt
(64)
-k.BB
-------
64
The concentrations within the slugs then become
A = A^At ,
(65)
From these expressions and Eq. (61), we find that the mean concentration
A is
A = A e-kAt = e-kAt t (66)
and a similar equation holds for B. Here, A,, denotes the initial value
of A~. Upon differentiating Eq. (66), we find that
- kA
dt ~KAA • (67)
Since this expression is not explicitly dependent upon £A or
the temporal behavior of A is independent of the way in which A is dis-
tributed within L. We can conclude, therefore, that linear reactions
are unaffected by the presence of concentration fluctuations.
Suppose, however, that A and B react with each other in the manner
A + B -* D
with rate constant k. In the example shown in Figure 20, it is clear
that this reaction can occur only in the domain (of width £Ag) where A
and B are mixed. Outside this domain, the concentrations of A and B
remain unchanged for all time. Let Am and Bm denote the concentrations
within the mixed zone £/\B- Then
dA dB
dT ' dT - -kAmBm ' (68)
-------
65
and from the definition given by Eq. (61), we have
A = U V£A - £AB> + VAB
Bm£AB
(69)
Keeping in mind that all terms on the right side of Eq. (69) are con-
stants except Am and Bm, we obtain upon differentiating Eq. (69)
= . . £AB
"
dt L dt "mm L ' (70)
and
= _kA B
dt KAm m
These equations show that in contrast to the linear reactions considered
earlier, the behavior of A and B is now explicitly dependent upon SL/\Q.
In other words, because of the nonlinearity of the chemistry, temporal
changes in A and B are sensitive to the manner in which the species of
A and B are distributed in space and time.
This fact is obvious from Figure 20, which shows that the quantities of
A and B that are consumed by the reaction A + B ->- D are determined by the ex-
tent to which the slugs of A and B overlap. Moreover, this figure indicates
that the extent of overlap, i.e. £/\B> is not expressible in terms of A and
F alone because these variables are independent of the positions of the
slugs within the interval L.
Herein lies the essence of the closure problem—the parameters that
control A" and F, in this case Am, Bm, and np^, are not expressible in
terms of only the variables A and B themselves. The origin of this diffi-
culty lies in the loss of information about the details of the distribu-
tion of A and B within L when A and B are averaged over the interval L.
-------
66
Consequently, in situations where these details play on explicit role
(namely, where nonlinear processes are involved), the closure problem
arises in equations governing mean values. This suggests that solving
the closure problem means retrieving the information that is lost in the
averaging process. This simple notion is the basic principle underlying
our proposed closure approximation. Specifically, restricting ourselves
to chemical reactions that do not perturb the fluid flow (all air polu-
tion reactions fall into -this group), we plan to show that most of the
information contained in the term A'B' is also present in its counter-
part A|Bj, which pertains to the case where the species A and B are
inert, and since AjBj is determined solely by those processes that con-
trol the spatial distributions of A and B, then AjBj, and hence A'B',
are expressible in terms of measurable properties of the system.
To demonstrate this approach, we consider the problem described in
Figure 20. Under the passive reaction assumption, the extents of the
regions 2,/\, SL$, and &AB are the same, regardless of the chemical proper-
ties of the species A and B. Consequently, if inert materials Aj and
BI are released into this system in initial concentration AQ and BQ,
and if the same mechanisms that controlled the positions and sizes of
the slugs of the reactive species A and B also act upon AI and BI, then
we would observe
— O^A
AT - —r» (71a)
and the second moments would be
-------
67
A2 VA .__ .
A = -— , (71c)
n
AIBI = L ' (71e)
All parameters on the right-hand side of Eqs. (71a) through (71e) are
identical to the corresponding terms entering into Eq. (69), and all
quantities on the left side of Eqs. (71a) through (71e) represent mea-
surable parameters. If the latter are regarded as properties of the
system and are specified along with the other parameters that describe
the given problem, then the five equations comprising the set can be
solved for the five unknowns Ag, Bg, £^, £R, and £^g. These results
can then be substituted into Eq. (69) to obtain the variables Am, Bm,
and £/\g that are required to achieve closure of the reaction rate equa-
tion [Eq. (70)].
The results of the solution of Eqs. (71a) through (71e) are
AB
J__
5
rA
(72)
*7
BI
R = —
-------
where
AB " A R
AIBI
(73a)
A
(73b)
BT
*"* J-
Fi
(73c)
After substituting these results into Eq. (69) and solving for Am and Bm,
we can write Eq. (70) in the closed form
dA
dt
krArB
LAB
A - A (1 -
AB
LAB
[A
(74)
Since the right side of this equation is equivalent to
-kAB = -k[AB + A'B1]
(75)
our approximation for A'B' is given by
A'B1 =
FArB
AB
A - AT 1 -
AB
B 'J
1 -
LAB
LA /J
- AB
(76)
We emphasize that this relationship is exact in the case of the problem
of two rectangular slugs of A and B shown in Figure 20.
-------
Before we extend this technique to more generalized situations, it
is of interest to test the accuracy of the Hilst et al. closure [Eqs.
(59f) through 59g)] in the context of this simple problem.
To simplify the mathematics, let us assume that
= a, =
„
= £,
69
i.e., that the species A and B are premixed. If A = B, then we find from
the definitions of the bar average ( ) and the fluctuation A1, Eqs. (61) and
(62), respectively, that
A
,3 _
+
+
(77)
For this quantity, the Hilst et al. approximation [Eqs. (59f) through
(59g)j gives
A'3
ADH
A3,
L
A3,
L
1 C ^ J
1.5 - L H
o c *• 4.
O - C> i T
2~
L2 J
2"
o ^
p
£ < 1/2
(78)
L/2 < a < L
Figure 21 shows the behavior of the fractional error of this approxi-
mation, i.e.,
e -
•5 -5
i J n i J
ADH
(79)
A
,3
as a function of the normalized mean square fluctuation
,
A2 "'
(80)
The large size of errors in this approximation, even for this elementary
problem, emphasizes the extreme difficulty of finding a closure scheme
based on hypothetical relationships among statistical moments that is
accurate for a wide range of situations.
-------
70
Consider now the extension of the closure scheme developed above
for the simple problem shown in Figure 20 to the generalized situation
exhibited in Figure 22. The disjointed volume v^g represents the total
volume where A and B particles are found together; and the volumes v^
and vn denote those regions where A and B are found in isolation. If
D
the averaging volume V implicit in the bar operator ( ) encompasses all
of Vn, Vg, and v^g, then we have by definition
I
A
A dv
LVA ~VAB
(81)
B = ^
B dv +
B dv
'AB
(82)
J
AB = -J- I AB dv
VAB
(83)
Let the following variables be defined:
i_ f
VA J
A - ~- \ A dv
= 77- B dv
(84)
a =
'AB
J A
'AB
dv
b =
'AB /
B dv
(85)
'AB
In terms of these variables, A and B become
-------
VI
oc
-o
LU
m o:
ii
u
FIGURE 21. FRACTIONAL ERROR e IN THE DONALDSON-HILST APPROXIMATION
[EQ. (59f)] OF A'3 FOR THE CASE OF TWO PREMIXED SLUGS OF MATERIAL IN
A ONL-PIMENSIONAL REACTOR (plotter! as 9 function of the
normalized mean square fluctuation
FIGURE 22. MIXING TURBULENT CLOUDS OF
REACTING PARTICLE SPECIES A AND B
-------
72
[AvA + avAB] , (86)
A = v|AvA+ avAB
= l|^BvB + BvAB] . (87)
Recalling that the concentrations within and the extents of the zones
VA and Vg are the same regardless of the chemical nature of the reacting
species, we find that the quantities corresponding to Eqs. (86) and (87)
for inert particle species ftj" and Bj are
Ai = Av + > (88)
BI
where
aT = — f AT dv , BT = — f BT dv . (90)
1 VAB J l l VAB J l
VAB VAB
Finally, we define the two parameters y and y such that
f AB dv = yabvAB , (91)
VAB
ATBT dv = yTaT6TvflR . (92)
I i 111 MD
VAB
-------
75
efforts to relate X to indices of both of the chemical natures of the
reacting species and of the turbulence and molecular diffusivity.
C. FORMULATION OF THE PARAMETERS ? AND r
M D
We consider three distinct situations, each of which illustrates to
some degree the physical significance of these parameters and how each
is related to properties of the flow.
1. Release of A into a Uniform Field of B
Suppose that a cloud of A particles of initial concentration AQ and
volume VQ is released into a turbulent fluid in which particles of species
B are uniformly distributed. We desire to estimate g\.
In this situation, v^, and hence CA, are initially zero. Although
particles at the edge of the cloud are near particles of the other species,
there is no region in which one species is immersed in the other. However,
with the passage of time, particles of A, excited by molecular-scale agita-
tions, migrate into the ambient fluid and thereby create a region v^g.
The flux of migrating particles at the cloud edge can be represented in
the form
(101)
where < is the molecular diffusivity of A and dA/dn is the gradient of
concentration of A particles normal to the cloud surface. If S denotes
the instantaneous surface area of the cloud, then at any instant the number
of A particles moving into the ambient fluid is given by
djl - f
dt ~ Jc
E • n ds (102)
S ~ ~
-------
76
Since the portion v^ of the original cloud that is occupied only by A
particles is completely surrounded by a region VAB created by the migra-
tion of A particles into the ambient fluid, all particles that leave v^
move immediately into v^g. Therefore,
- N
where NQ = AnVg is the total number of A particles released and N is the
number still in v^ at any given time. From Eqs. (103) and (102), we
obtain
= F • n ds (104>
dt F S
The particle flux F averaged over the cloud surface can be expressed
in the form
F = < - , (105)
where A and £ are concentration and length scales, respectively. We pro-
pose that the ratio A/? entering in Eq. (105) can be represented to good
approximation using the following expressions for A and 5 :
- , (106)
V0
•$
where e is the energy dissipation rate of the turbulence.
-------
77
The quantity on the right side of Eq. (107) represents the smallest
length scale associated with spatial variations in a scalar field in
turbulent fluid (see Batchelor, 1967). It represents the length scale
associated with the concentration gradient that is required to produce a
molecular diffusive flux of particles that just balances the compressional
fluxes arising from the strain rate field of the turbulence. This length
scale is thus a logical choice for £, at least initially. With the passage
of time, however, ? increases, but in a manner that is not possible to
describe exactly. We are proposing that the overall effect of this in-
crease (and the simultaneous decrease in A) on the flux F can be accounted
for by Eqs. (105) through (107); i.e.,
- . (108)
0
Here, VQ is constant because it is the combined volume occupied by all
particles that initially were of species A. Combining Eqs. (108), (103),
and (104) , we obtain
dcfl
It is not possible to obtain a precise description of the rate of
change of surface area S of a cloud of marked particles in turbulent fluid.
A reasonable approximation for the average surface area might be
S ^ a2 , (110)
where a is the distance between a pair of marked particles released a
distance £Q apart where
-------
78
and Sg is the initial surface area of the cloud. Batchelor (1950) showed
that for particles whose separation I is within the inertial subrange of
a homogeneous turbulence,
d£2
Combining this result with Eqs. (110) and (111), we obtain the estimate
S = S0(l + £2/3 SQ-2/3 t2) ; (112)
and upon substituting this expression into Eq. (109) and solving the re
sulting differential equation, we obtain finally
, (113)
where
a = 3K1/4 e'/4 ^ , (114)
b - .'/" c11/12 *-7/3 . (115)
and £Q is the initial cloud diameter.
Under conditions of strong convection in the atmosphere with e on
the order of 10~2 m^/sec3, Eq. (113) predicts that for a cloud of material
of Schmidt number on the order of 1 and initial diamter of roughly 1 meter,
CA will reach its limiting value of unity within about 20 seconds after release.
-------
79
2. Release of Dilute Mixtures of A and B
into a Homogeneous Turbulent Reactor
In view of its definition [Eq. (96)], c^ is equivalent, in a homo-
geneous system, to the probability that any inert particle of A will be
found at time t within some distance 3r, say, of a particle of B. In other
words, implicit in our definition of v/\g is the condition that the smallest
distance between particles of opposite species is on the average 3r or
smaller. The magnitude of 3r is in turn an implicit function of the chemi-
cal reactivity of A and B because reactions are assumed to occur within
and only within
Let #|AQ represent the event
E |A« = given particle A~ is within a distance
8r of at least one B at time t . (116)
If the concentrations of A and B are sufficiently dilute that the proba-
bility of finding more than one B particle with a distance 8r of the given
A is negligible, then clearly
prob B|A0] = //// P(r»t,r',t|rAO,tAO,rBO,tBO)SB(rBO,tBO) dr' dr dtBQ dr
0 0 8V
(117)
where (r/\o»t^Q) is the release point of the given inert A particle, 8v is a
volume of radius 8r centered at the location (r,t) of a particular inert B
particle, and SB(rBQ,tgg) is the number density -of B particles released per
unit time and volume at
Equation (117) gives the probability that a particular A particle ex-
periences event E. However, CA is equivalent to the event
-------
80
E = any randomly chosen particle of A is within a distance
6r of at least one B at time t.
If N/\ particles of A have been released prior to time t, then there are
N/\ mutually exclusive ways in which E can be realized, and it follows that
il i I f \
= H~ II Pr0'") ^l^n ^A^rAO'^AO^ ^AD ^rAO ' ^^^
A 0
where S^(r^Q,t^Q) is the number density of A particles released per time and
volume at (rAO^o)- Upon substituting Eq. (117) into Eq. (118), we obtain
finally (see Lamb, 1974)
AT(r',t) BT(r,t) dr1 dr . (119)
XA "^ ... L ~ i ~ ~ ~
If the volume 6v is sufficiently small that Aj(r',t) Bj(r,t) is nearly
constant for all r' in 6v and for all r, then Eq. (119) reduces to
Aj(r,t) Bj(r,-
dr . (120)
A
If, moreover, the chemical reactor has volume V and the statistics of Aj
and BI are uniform throughout V--the homogeneity assumption—then Eq. (120)
gives finally
. (121)
-------
81
In the limit as A and B became thoroughly mixed, r^g ->• 1 and ^ ^ BI 5v<
Since by definition 0 < r,A <_ 1, we see that Eq. (121) holds only if the
concentrations are so dilute that
FT 6v < 1 (122)
3. Release of Rectangular Slugs of Reactants
into a Hypothetical, One-Dimensional Reactor
Imagine a one-dimensional system such as that shown in Figure 20 in
which the reacting clouds have rectangular shapes. In this hypothetical
system, we see that
A0£AB _ £AB _ £AB
where AQ is the (uniform) concentration o^ inert particles in the cloud.
The following mean values are thus defined:
A-
A
I L
A£n R^j)
.2 VA 2 _ VB
Solving these six equations for £y\g and a^ and substituting the result into
Eq. (123), we obtain
-------
82
and, similarly,
where
(125)
= —
A -2
AI
- -
B"-2
BI
AB ATBT
It turns out (as we show later) that when Eqs. (124) and (125) are used
in the closure approximation [Eq. (95)], good agreement with observation
results, even in turbulent pipe flow reactors. This suggests that Eqs. (124)
and (125) may be sufficiently accurate for all general purpose applications
of Eq. (125), but this remains to be definitely established.
D. FORMULATION OF THE PARAMETER A
Being dependent upon the properties of both the turbulence and the
chemical reaction, A is more difficult to formulate than the mixing para-
meters C ar)d C- By virtue of its definition, i.e.,
-------
83
AB dv
1
'r /A
dv
AB
B dv
AB
Bj dv
(126)
IjBj dv
AB
X is a measure of the effect of the chemical reaction on the correlation
of the concentrations of A and B within the mixed zone V/\B- This can be
seen more clearly in Figure 23, which compares profiles of the concentra-
tions of inert species in v/\g with those that might be observed if the
species were reactive. In interpreting this diagram, one should keep in
mind that v/\g is defined for the inert rather than the reactive particles.
As a consequence, portions of the zone defined as v/\g may contain only one
of the species rather than a mixture of the two. In fact, if the reaction
is extremely fast, no part of V/\B contains a mixture of the reactants ex-
cept a very small zone where the reactant clouds merge. In this case, it
is clear from Eq. (126) that A -> 0. However, at the other extreme where
the reaction is exceedingly slow, the concentration profiles within v/^
are nearly identical for both the reactive and inert materials, and it follows
that A -»• 1.
- X
FIGURE 23. COMPARISON OF CONCENTRATION PROFILES IN THE MIXED
ZONE VAB FOR INERT AND REACTIVE SPECIES
-------
84
The behavior of the parameter x for reactions of intermediate speed
has been studied by Shu (1975). In his analyses, the differential equa-
tions governing the combined reaction and molecular diffusion of two reacting
species in one dimension were solved numerically, and the values of X were
calculated explicitly for a wide range of rate constants, molecular diffu-
sivities, initial conditions, and the like. The model equations used were
9A
at
/A
- kAB
(127)
9t
A
- knAB
(128)
(129)
8t
(130)
where k is the reaction rate constant, D is the molecular diffusivity, and
n is the stoichometric ratio. These equations were solved subject to the
initial and boundary conditions
A(x,0) = Aj(x,0) = <
AO , -6 < x <_ 0 ;
0 , otherwise;
(131)
B(x,0) = Bj(x,0) = <
, 0 < x <_ 6 ;
0 , otherwise;
(132)
lim A, B = 0
(133)
-------
85
The computed functional form of x(t) is shown graphically in Figure 24 for
several values of the dimensionless group
tD k62nA
a = ^= ——- , (134)
LR u
where tD is the characteristic time scale of the diffusion; i.e.,
(135)
and tg is the (initial) characteristic time scale of the reaction,
namely,
tR = (knAor . (136)
All curves in this figure are for the case of the feed ratio
(137)
equal to unity.
It can be seen that in all cases both the initial value and the final
values of A are 1. Initially, before the reaction has had time to consume
appreciable quantities of the reactants, the concentration profiles are
virtually identical to those in the inert case. In the long time limit, the
eventual total mixing of the species results in uniform spatial concentra-
tion distributions that, as is easily verified, are characterized by x = 1,
regardless of the value of k.
-------
It is important to note that for a < 1, the parameter x has approxi-
mately unit value for all time (see Figure 24). This greatly simplifies
the functional form of the closure scheme [Eq. (95)]. In contrast, signi-
ficant temporal variations occur in X as a becomes large compared with
unity. A fortuitous and fortunate feature of this behavior is that for all
values of a, X reaches its minimum value Xm-jn at approximately t = 0.5tp and
maintains this value for a period of about 100 tp. Inasmuch as reacting
species that are initially unmixed can be consumed only as fast as the mole-
cular diffusion and turbulence can bring them together, for fast reactions
(a » 1), the characteristic time scale must be tp rather than t^. In view
of this, it appears from the behavior of x(t) described above (and shown in
Figure 24) that a constant value of X, namely,
X - »mjn , (138)
should yield good results in Eq. (95) for all values of a. With this pros-
pect in mind, Shu examined the relationship between Xm-jn and a and found
that for 6=1,
1 T T2"
as shown graphically in Figure 25.
An interesting consequence of Eqs. (138) and (139) is that for large
values of a, we have
Upon substituting this result into Eq. (95) and using the resulting ex
pression for JB in the rate equation for A, we obtain
-------
10
-1
10
-2
10
,-3
a = 1
To"
100
1,000
10,000
10
-4
Source: Shu (1975).
1.0"
10"^ 10"'
t/tD = t/(«tR)
10
0
10'
FIGURE 24. FUNCTIONAL FORM OF X(t) FOR VARIOUS VALUES OF THE RATIO a
OF THE DIFFUSIVE AND REACTIVE TIME SCALES tD AND tR
oo
-------
10
•icr1
10
-2
10
-3
10
-4
10
-1
10'
Source: Shu (1375).
10'
a
FIGURE 25. DEPENDENCE OF X . ON ex FOR THE CASE OF THE FEED RATIO (3 = 1
CO
en
-------
89
0
The absence of the rate constant k from the right side of this equa-
tion emphasizes the fact that fast reactions are diffusion-limited processes.
Shu (1975) also considered the effects of variations in the feed ratio
3 and the number and initial separation of the slugs of reactants. Figure 26
shows how x(t) is affected by variations in 3 when a is held constant at 700.
The most striking effect is the more rapid onset of the rise of x toward its
final value of 1. This behavior reflects the more rapid exhuastion of the
species (in these cases A) that is present in a smaller quantity. Despite
this behavior, the assumption
x = constant = A .
should still prove accurate during the major part of the period in which
the reaction is in progress. Figure 26 also shows that except for the case
of 3 « 1, the order of magnitude of Am-jn is unaffected by changes in 3.
Figure 27 shows the influence on x of changes in the number and sep-
aration of the initial slugs of reactants. In each of the three separate
cases shown, the minimum value achieved by X is of the same order of magni-
tude. These results, those presented above, and other experiments not des-
cribed here suggest that Eqs. (138) and (139) are perhaps accurate even in
general situations. In the next section, we examine this possibility in
light of observational data.
E. TESTING THE CLOSURE SCHEME REPRESENTED BY EQ. (95)
USING OBSERVATIONAL DATA
Few empirical studies of turbulent reactions have been reported in"
which all of the information required to apply and subsequently to test
Eq. (95) has been measured. Usually, the data reported are not sufficient
-------
10'
Source: Shu (1975)
FIGURE 26. IMPACT OF VARIATIONS IN THE FEED RATIO B
ON THE FUNCTIONAL FORM OF x(t) FOR a = 700
-------
10° H
iif'h
10
-2
10
-3
10
-4
10
-5
10
-4
10"
Source: Shu (1975).
10
-2
10
-1
T_
Tn
10
10'
10-
FIGURE 27. INFLUENCE OF CHANGES IN THE INITIAL NUMBER AND
SEPARATION OF REACTANT SLUGS ON THE FUNCTIONAL FORM OF X
-------
92
to permit estimations of the parameters FAB> CA» and <;g that are properties
of the inert species concentration statistics. However, using the theory
of Toor (1962), which relates the concentration statistics of inert species
to those of extremely fast reactants in the same turbulence, Shu (1975)
succeeded in collecting experimental data from the literature that are ade-
quate to apply and test Eq. (95). (See Appendix E for details.) We report
some of his findings here. For more details, we refer the reader to Shu (1975)
All of the data used were gathered from the multijet reactor designed
by Vassilatos and Toor (1965). In this system, the reactants, usually an
acid and a base, are fed through alternate jets in a head consisting of
100 small jets at the end of a pipe several centimeters in diameter. Down-
stream from the jet head, the mean concentrations are constant on all planes
normal to the pipe axis. Thus, this reactor approximates a one-dimensional
system in which downstream distance corresponds to elapsed time.
Using Eq. (95) in the governing rate equations, Shu modeled this
system by
dFA kArAR__
— - fl R TF C\ - r \~\ TF n r M ("\ A.~\ }
P TT ATBT iT-n ~ U t>i\)l irD - U - CrJ J , U^U
dt rnrn nlul L'A v' -°A'J l'B
where
FA = , (142a)
AI
FR = , (142b)
D n"
BI
and where Aj and Bj are constants. Since the reactions are of the form
A + nB -»• products, and since the mean speed IT of material down the pipe
is constant, Eq. (141) can be written in the equivalent form
-------
93
dFfl k*r nAT
3T = --^§~1 [FA - (1 - ?A>] [FA - ] + ^ ' 043)
where x represents downstream distance and
is the feed ratio of the inlet jets.
Shu used Eqs. (124) and (125) to approximate ^ and eg. Subsequently,
he determined the values of f^, fg, and f^g entering into these expressions
and into Eq. (143) by using Toor's theory, mentioned above, together with
relevant measured data on extremely fast reactions in the multijet reactor
(see Appendix E). We should add that although Toor's theory has been corro-
borated by several experiments, a more thorough and meaningful test could be
performed if data were available that did not require this theory.
To apply the approximation [Eqs. (136) and (139)] for A, Shu evaluated
a using reported values of k, D, and other terms, and the measured scalar
microscale for 6. In all cases, Shu found that a £ 1, and so he used the
constant value x = 1 throughout the model validation calculations.
Figure 28 presents the first set of comparisons between the observed
values of F^ and those predicted by the model equation [Eq. (143)] for two
different values of the feed ratio B. In both cases, the agreement between
the model predictions (solid line) and the observations (circles) is very
good. The degree to which the concentration fluctuations dominate the
rate of change of the mean reactant concentrations is revealed by the dis-
crepancy between the data points and the dashed lines shown in both panels
of Figure 28. These lines show the predictions of the rate equations in
which the concentration fluctuation term TrE"1" has been omitted. The com-
parisons show that the effect of the fluctuations is a retardation of the
-------
CO
t
o
CO
o
X
O
OJ
rr
o DflTn OF MflO. RE=2U30
- D. Z. MODEL
k = 0.83 x TO4 I mole"1 s"1
o
o
FLUCTUATION
OMITTED
o
I I
0.0
1 .0
2.0
X,CM
Source: Shu (1975).
3.0
.o-
0.0
j j } r
o ORTfl OF MRO. RE=2H30;'
- D. Z. MODEL
6 = 1.5
U.O
FIGURE 28. FIRST COMPARISON OF CLOSURE MODEL PREDICTIONS WITH OBSERVATIONAL DATA
-------
95
speed of the reaction. This result is caused by the segregation of the
reactants initially and the finite time required for the turbulence and
molecular diffusion acting together to bring the two materials into contact.
As shown in Figure 29, this phenomenon cannot be described in terms of A and
F alone. In this figure, we have plotted the solutions of the equation
rtT .
31= -ekAB (144)
for several values of the constant e. Note that no value of e exists
that will bring the solution of Eq. (144) into even crude agreement with
the observations.
Figure 30 shows a comparison of the theory and observations for a situa-
tion involving both a faster reaction and larger values of the feed ratio
than those presented in Figure 28. As this figure indicates, the agreement
between theory and observation is better than it was in the previous case.
Additional comparisons can be found in Shu (1975).
The results reported above are very encouraging and provide the hope
that our closure model will perform well under a wide range of conditions,
both in the atmosphere and in the laboratory. Our immediate goal is to
develop expressions for r/\g, c/\> and Cg that permit Eq. (95) to be applied
to reaction processes in the atmosphere. Toward this end, we are currently
preparing a numerical experiment in which these functions will be estimated
using the turbulence model of Deardorff. The procedure will be to release
pairs of particles in the numerical fluid in each of many realizations of
the turbulent flow and to calculate the mean square concentration field from
the ensemble of particle-pair trajectories generated thereby. In this
way, we plan to derive the forms of f^g, c/\> and eg as functions of the
atmospheric stability, wind speed, release height, surface roughness, and
the like. Provided that the resulting expressions have universal features,
the universal functional forms of these parameters can be used in conjunction
with those of the diffusivities KZ presented in Chapter II to study reactive
plumes of material in the planetary boundary layer. Progress on this project
will be reported in Shu (1975) and in future reports on our EPA contract
efforts.
-------
7.0
A OBSERVATIONS
k * 8320
0 - 1.5
0
FIGURE 29. ATTEMPTS TO SIMULATE THE OBSERVATIONAL DATA USING A RATE EQUATION
THAT OMITS CONCENTRATION FLUCTUATION EFFECTS
CTl
-------
VflSSILRTOSClSSU)
K = 0.12U.E 05'
BETfl = 1.Z6
0.0
4.0
o
CO
o
CO
o
X
CM
O
O
O
0.0
.0
VRSSILRTOSC19GU)
K = 0.12UE 05
BETfl =3.98
2.0
X,CM
3.0
4.0
Source: Shu (1975),
FIGURE 30. SECOND COMPARISON OF CLOSURE MODEL PREDICTIONS WITH OBSERVATIONAL DATA
-------
98
IV DEVELOPMENT OF A SCHEME FOR PARAMETERIZING THE EFFECT
OF SUBGRID-SCALE CONCENTRATION VARIATIONS
ON REACTION RATE
A. INTRODUCTION
The previous chapter dealt with the instantaneous concentration fluctu-
ations produced by turbulence and the effect they have on the reaction rate
of nonlinear chemical reactions. In this chapter, we treat the analogous
problem of variations in the mean concentration that are of too small a scale
to be resolvable by a grid model of an urban atmosphere. The similarity of
this problem to that of the turbulent concentration fluctuation is so great
that the closure scheme we developed in the last chapter can be used here as
the basis of a parameterization scheme for the subgrid-scale effects.
Before proceeding with the development and application of the scheme,
we briefly review the nature of the subgrid-scale variation (SSV) problem
and then develop necessary conditions under which SSV is of negligible con-
sequence. The last two sections of this chapter are devoted to the deriva-
tion of the SSV parameterization scheme and to the application of it to air
pollution problems of applied interest.
B. STATEMENT OF THE PROBLEM
Consider two pollutants, A and B, that undergo the reaction
A + B -> D
to form a third product D. The equations governing the concentrations of
these species in an urban atmosphere are generally taken to be (for modeling
purposes)
-------
99
where
A = some time averaged concentration of A,
u- = the i-th component of the mean wind speed,
K. . = the turbulent diffusivity tensor,
' J
S,, = the distribution and strength of sources of A,
k = the rate constant for the above reaction.
An equation similar to Eq. (145) governs the concentration of B:
^D f\D r\ ^D
(145b)
Because the analytic solution of this system of equations is not available,
Eq. (145) must be solved numerically. For reasons given in Chapter I, this
necessitates the averaging of Eq. (145) over volumes equal to that of the
grid cells into which space is discretized in the numerical technique. The
spatial average is denoted by the tilde (~) and is defined by (using the con
centration of A as an example)
X+AX y+Ay Z+AZ
x L L
where r = (x,y,z) and (ZAX, 2Ay, 2Az) denotes the dimension of the grid cells.
Averaging Eqs. (145a) and (145b) in the manner of Eq. (146) and assuming that
u and K are nearly constant over distances comparable to the grid cell dimen-
sions, we obtain
where SA and AB are defined as in Eq. (146). A similar equation can be derived
for B. Equation (147) and the corresponding equation for B are used as the
basis of urban pollution models; therefore, it is A and B, rather than the
"point" values A and B that these models predict. The differences
-------
TOO
A" = A - A
(148)
B" = B - B
are the subgrid-scale variations. Using the definitions given by Eq. (148)
and assuming that
A = A
B = B
and so forth, we can express Eq. (147), and the corresponding equation for
B, in the more useful forms
||+ LA = SA - kAB - kA^B" , (149a)
|| + LB = SB - kAB - k^B" , (149b)
where for brevity we have introduced the linear operator
L -= =1 157 - it "u it • "50>
I I J
The quantity kA"B" on the right side of Eqs. (149a) and (149b) represents
the influence of the SSV on the rate of decay of A and B. Until now, all
models of urban pollution have tacitly ignored SSV effects, but as we show in
the next section, this omission is generally unjustifiable.
C. CONDITIONS UNDER WHICH SSV EFFECTS
ON REACTION SPEED ARE NEGLIGIBLE
If
A"B" « AB , .. (151)
then it is clear from Eq. (149) that SSV has a negligible influence on A and B
and can be ignored. To derive the condition under which Eq. (151) is satis-
fiedv we begin by deriving the equations governing A" and B".
-------
101
Subtracting Eq. (149a) from Eq. (145a) and using Eq. (148) to define
A", we obtain
j.nll - _ , J
— ^ LA" = SA' - k(AB" + BA" + A"B" - A"B") , (152a)
and, similarly,
P.DII ~ -, , <•
H + LB" = Sg - k(AB" + BA" + A"B" - A"B") , (152b)
where SA' = SA - SA and S£ = SB - Sg. Multiplying Eq. (152a) by B" and
Eq. (152b) by A" and adding, we obtain
3 fA"R"^ 4. rfA"nM - v I 9 A 3B , 3A 9B
-(A B ) + L(A B ) - - K. . — -^- + ^- ^^
L- ' J J I _
SA'
- k(AB"2 + BA"B" + A"B"2 - B"A"B"
+ AB"A" + BA"2 + A"2B" - A"^")
Upon averaging this equation in the manner of Eq. (146) and invoking assump-
tions such as A = A made earlier, we obtain
. r^
'ijL8xi
ill ~AII ~nii
Annil I rn'iDli - V I on OL> ' °" °t>
B ~t~ I/A 15 - - l\ . .
1J
J J U
- k[AB"2 + BA"2 + /TB" (A + B) + A"B"2+ A"2B"]
(153)
The first term on the right side of this equation causes a decay of A"B",
even in the absence of chemical decay, i.e., when k = 0. Thus, if A and B
were inert and if A"B" were nearly uniform in space, the steady-state value
of A"B" would be such that
L3_A" 3B" , 9A" 8B" - RHC" j. /\'ic»
fe--37: + ^r:^7 -B SA + ASB
L. ' J J ' -J
-------
102
r ^^~ "^
Let the steady-state value prescribed by Eq. (154) be denoted by (A"B")T.
It follows, therefore, that the conditions that result in the satisfaction
of
(f-B-'Jj « AB (155)
are sufficient conditions for the realization of Eq. (151). Since these
conditions are not difficult to derive, we use them as the measure that we
seek of the importance of the SSV.
We note first that the left side of Eq. (155) can be expressed approxi-
mately in terms of known quantities using Eq. (154) and the equation corre-
f^~~i*>
spending to Eq. (153) that governs A" , namely,
aA" ? aA" aA"
on + rA" = -9Y. + A"c;" - KAA"R" + RA" 4- A"
3t LA ^ 1J3X. 9X. A K(AA B BA A
(156)
In Appendix F, we show that the first term on the right side of Eq. (156)
can be approximated by
\ ' 9
A • <157'
where X is on the order of the smallest dimension of a grid cell and K^ is
the turbulent diffusivity in the corresponding coordinate direction. Also,
it is not unreasonable to assume that A"Sj! is on the order of
(158)
From this approximation and from Eq. (157), we conclude from Eq. (156) that
in the steady state
£ , (159)
-------
103
ignoring the chemical decay (k = 0) as before. Employing assumptions such
as Eqs. (158) and (157) in Eq. (154), we conclude finally that
The quantities on the right side of Eq. (155) pertain to the concentra-
tions of the reactive species A and B, rather than to their inert counterparts
considered above in deriving Eq. (160). Since the decay rates of A and B are
on the order of
-TV- = -rr '^ -kAB , (161)
and since the production rates of these quantities are SA and SB, respectively.
we see that in the steady state
AB * T^ (a? + ?)' , (162)
K B + A
where
a = ^ . (163)
Assuming that
we can reduce Eq. (162) to the form
2SflSR
AB A B
k(sA + SB)
Upon combining this result with Eq. (160), we conclude finally that a
sufficient condition for ignoring SSV effects on the chemical decay of A
and B is
-------
104
2K ?ASB
~
k(S
(164)
When this condition is not satisfied, SSV effects may or may not be important,
depending on the particular physical situation. In the next section, we con-
sider the derivation of a necessary condition under which subgrid-scale
concentration variations significantly influence A and B.
To demonstrate the quantitative significance of Eq. (164), let us consider
a simple one-dimensional problem in which a series of sources of strength S and
width W are distributed with separation £ along the x-axis, as shown in Figure
31. We might consider this problem to represent a cross section of a system of
infinitely long, parallel streets.
T
S
I
FIGURE 31
HH
SYSTEM OF RECTANGULAR SOURCES
Suppose that S represents the emission rate of species B and that aS is the
emission rate of species A. Suppose further that a grid model simulation of
this problem is to be constructed using a grid mesh size Ax. We desire to
know whether subgrid-scale concentration variations can be ignored in the grid
model simulation of the concentrations A and B.
Through a straightforward series of analyses, we find that in this problem
S"S" =
VB
and that
I
1 -
~
S =
aSw
I
(165)
(166)
-------
105
Substituting these results into Eq. (164), we find that SSV can be
ignored if
2K?
kS < -, - * - - . (167)
Normally £»w and y might be on the order of one. Under these conditions,
Eq. (167) can be simplified to
<
kS $ ~- . (168a)
A
For the hydrocarbon-ozone reaction, k ^ 0 .1 ppnf min" . Also, hydrocarbon
emissions from motor vehicles are on the order of 10 gin/mile. Using these
values we find that for a roadway carrying 1000 vehicles per hour,
kS * 3 x 10"5sec"2 . (168b)
2
A typical value of K in the atmosphere is 10 m /sec. Consequently, we find
A
from Eq. (168a) that SSV effects are negligible in the present problem if
A ~ 50 m . (169)
In the problem shown in Figure 31, A is on the order of £. Thus, if the
source configuration shown here represented a network of actual streets, £
would certainly exceed 50 m, and SSV effects could not be ignored. In cases
where £>Ax, A is on the order of Ax. In this case, SSV effects could be sup-
pressed by making the mesh size Ax suitably small; but this is not practical
in many cases, and one is left with a problem in which SSV effects might have
to be parameterized to achieve meaningful estimates of the spatial averaged
concentrations A and B. In Section E, we develop a quantitative measure of
whether parameterization is actually necessary.
-------
106
D. PARAMETERIZATION SCHEME FOR THE SUBGRID-SCALE
CONCENTRATION VARIATIONS
Having derived the conditions under which SSV effects are negligible
in the modeling of nonlinear chemical reactions and having demonstrated
that these conditions are frequently not realized in actual air pollution
simulation studies, we turn now to the development of a mathematical scheme
for representing SSV. The scheme we developed here is based on the closure
approximation developed in Chapter III.
Let V denote the volume of any grid cell in the grid model of A and B,
and let v,, denote the total volume within V in which particles of species A
are found alone. Similarly, let VB denote the volume where B particles are
found alone. The volumes vfl and VR are distinguished from a third volume
M D
VAR in which particles of A and B are mixed sufficiently for chemical reac-
tions to occur. Figure 32 illustrates these three volumes. By definition,
we thus have
A dv +
A dv
(170a)
'AB
/ r c
B = 1 [J B dv + J
B dv,
(170b)
'AB
FIGURE 32. DESCRIPTION OF THE VOLUMES v., VB> AND VAB
-------
107
It is clear from the definition of v* and vg that if A and B were chemically
inert, the concentrations within these volumes would be the same as in the
reactive case. It is only in the domain VAB that the inert and chemically
reactive systems differ. This fact is one of the key points upon which our
parameterization scheme is based.
We define the following four variables:
A EE 1 /Adv , B = - /
VA J VB J
VA VB
B dv
a = 1 /A dv , b EE 1 /B dv . (171)
VAB J VAB J
VAB VAB
Let the subscript I denote the situation where the given particle ?pecies
is chemically inert. Then analogous to Eq. (171), we define
^T = TT f AT dv » h ET7 /" BT dv • (172)
1 VAB J l VAB y
VAB VAB
In terms of these variables, A, B, Aj, and B, become
(173a)
(173c)
AB] ' ' "073d)
-------
108
Treating v^, Vn, and VAB as invariants in the inert and reactive systems
as we have done above is permissible as long as the fluid motions, which con-
trol these volumes, are not perturbed by the chemical reactions. In the
problems of interest to us, this condition is always satisfied.
Next, we define
'AB
(174a)
A2 = 1
AI V
A2 dv + f A2 dv
v
AB
+ V!VAB]
(174b)
r =
rB -
B2 dv +
J B2 dv
VAB
.if
M
yRB2v
D
(174c)
(We discuss the p's later.) Just as Aj and Bj are determined completely
by the fluid velocity field and source distribution, so too are r^g, r/\,
and TQ. Consequently, these five quantities can be regarded as known
parameters that characterize the given problem.
Rather than parameterize the term A"B" explicitly, we shall instead treat
the term
AB = A"B" + AB
(175)
Using the variables introduced above, we can write
AB = 77
AB dv =
(176)
'AB
where x is a parameter, that can vary with time. The desired parameterization
is then obtained by examining a, b, and vnc in terms of the seven known quan-
Ab.
tities given by Eqs. (173), (174), and by substituting the results into Eq.
(176). To accomplish this, we must use two additional approximations because
-------
109
there are two more variables in Eqs. (173) and (174) than there are equa-
tions. We therefore write
aj = ?AA (177a)
bz = ?BB . (177b)
In terms of the new parameters 5. and cR5 AT and BT become
-------
no
where
A ~ ~~0 ^D = ^o~ r /,D - —r^r . (181)
A A2 , B B2 , AB AB
Using all of the above expressions, we obtain finally
(182)
where
h = -*- . (183)
Except for the parameters h and g, which we discuss in more detail later, all
of the terms on the right side of Eq. (182) are known quantities. The next
step, therefore, in implementing this parameterization of AB is to express the
right side of Eq. (182) in terms of the source distribution, fluid velocity,
and other variables that characterize the given system.
This task is simplified somewhat in situations where the species A and B
are released from the same sources, as is often the case in actual air pollu-
tion simulation studies. Under this condition, we have
'A = VB
and gA and gD reduce to
H D
11 11 „ _
(184)
Using these results, we find that
s* ^ ^ \
(185)
-------
m
A similar expression can be found for i\g,, - r.g. Now, from Eqs. (177) and
(178), we conclude that the term in parentheses in Eq. (185) is identically
zero. Thus,
AD
which reduces with the aid of Eqs. (180c), (174a), (178a), and (178b) to the
final form
f*^f ^ n, ^
AB - hrABAB . (186)
This simple form of the parameterization is valid whenever the reacting species
are emitted from the same source or in other terms, when the reactants are pre-
mixed. That most of the air pollution modeling problems of interest fall within
this class of problems is justification for our continuing the analysis of the
parameterization scheme with Eq. (186) rather than the general expression Eq.
(182). Thus, in the next section, we consider h briefly and then derive ex-
pression relating f«R to measurable parameters.
E. DERIVATION OF ?AB IN TERMS OF MEASURABLE PARAMETERS AND
APPLICATION OF THE PARAMETERIZATION SCHEME TO SPECIFIC PROBLEMS
Before considering i\B, we should point out that h is also a function of
the given system and ranges from a value of zero for infinitely fast reactions
to a value of unity for inert materials. Chapter III further analyzes this
parameter. In the present analysis, we assume that h is a known constant.
By virtue of its definition, fAB pertains to- the inert species Aj and Bj
and can accordingly be derived from linear theory. That is,
AB : -
Ari
-------
112
where AT is the solution of Eq. (145a) with k = 0; i. e. ,
9AT
-+ £A = S , (187a)
and, similarly,
3BT
+iBS . (187b)
We assume that A and B are released from the same sources. Hence,
SA = aSB , (188)
where a is the ratio of the emission rates of. A and B. If the initial con-
centration of AT and BT are both zero, then it follows from Eqs. (187) and
(188) that
and hence that
/V> ~
2 2
fAR = 4--=-4 . (189)
R A
BI AI
The solution of Eq. (187a) is of the form
Aj(r.t) = jfj p(r,t|r',t') S(r',t') df dr' , (190)
where p is the Green's function of Eq. (187a) and where, for brevity, we
have dropped the subscript from the source function S«. All aspects of
the transport, diffusion, and interactions of the particles with bound-
aries are embodied in the kernel p. If the turbulence is homogeneous,
then p has the property
P(r,t|r',t') = p(r-r',t,t1) , (191)
-------
113
which leads to the following simplifications of the expressions for AT and
""?
.t
p(rrr',t,t') S(r',t') dt' dr' dr, , (192)
v(r) u
ff p(r-r',t,t') ?(r',t') dt1 dr' , (193)
t.t
)(r^r',t,t') p(r-r",t,t") S(r',t') S(r",t") dt1 dt" dr1 dr" dr],
(194)
~ rr /*t t
Aj(r,t) =11 I C p(r-r',t,t') p(r-r",t,t") 6(r',r",t',t") dt' dt" dr1 dr"
JJ JQ JQ
(195)
where
>.AX ,.,Ay_ ~Az
6(r'fr",t'.t") =1| 2 12 JT S(r'+c,t') S(r"+s,t") d5
~ ~ Ifi ~~ ~~
J-AX_ J-AV_J-AZ_
222
(196)
and V = AxAyAz is the volume of the grid cell v(r) centered at r. The derivations
of Eqs. (193) and (195) are presented in Appendix B. Since p is the known Green's
function of Eq. (187), the parameter r\n can be evaluated using Eqs. (193) and
(195) once the source correlation function G has been prescribed. (Note that S
is the source function that enters into the grid model of A and B.) We consider
below the forms that G acquires in three problems of practical interest.
-------
114
1. Random Distribution of Point Sources in a
Three-Dimensional Space
Consider the situation in which "Point" sources of (small) volume v are
distributed at random in space, and suppose that each source has a constant
emission rate m. Let the number density of sources D(r) be sufficiently
large that on scales comparable to that of the grid cells used in a numeri-
cal calculation, the density D can be treated as a continuous function. The
problem just posed might serve as a model of emissions arising from home
heating units, or, through proper choices of D, from clusters of major point
sources of pollution.
The source strength function that describes this array of sources is
•
S(r,t) = ™ U(r) , (197a)
where
!1 , if r lies within any source ;
(197b)
0 , otherwise
Thus,
z+Az y+Ay x+Ax
S(r) = 1 f f J S(r') dx' dy' dz'
Z-AZ ^y-Ay x-Ax
= mD(r) , 098);
where we have taken the averaging Volume V to have the dimensions 2Az • 2Ax • 2Ay
Similarly, we have
AZ AX Ay >2
G(r'>r")=}f I I \U(r' + 0 U(rn + ?) dsy d?x dsz . (199)
•'-A 7 -AX -AV
-------
115
Let
W(r,Ar) = U(r) U(r + Ar) . (200)
It can be seen from the definition of U [Eq. (197b)] that W is a random
variable whose value is either 1 or 0 depending on whether both r and
r + Ar lie within sources. If we let
Ar = r1 - r" , (201)
Eq. (199) then can be written in the form
G(r',r") =
m
v2
X"+AX y"+Ay Z"+AZ
J J J
•'v/" Aw «'i<" Aw */-*
W(r,Ar) dr
x"-Ax y -Ay z"-Az
(202a)
or
• 2
G(r',r") = ^W(r",Ar) , (202b)
v
where W(r",Ar) is the Volume average of W at the point r".
Provided that the source density D in the vicinity of r" is sufficiently
large, we can exploit the random variable nature of W described earlier to
write
W(r", Ar) = (1) P(l;r") + (0) P(0; r")
=P(l;r") , (203)
where P(n;r") is the probability that W(r",Ar) has the value n, where n = 0, 1.
If Ar = (Ar . Ar , Ar ) is larger than the dimensions of the elemental sources,
v x y z 3
i.e.,
Arx > vx •
Ary > vy
(204)
-------
116
where v = v v v , then P(l;r") is just the product of the probabilities that
x y i.
the vectors r" and r" + Ar fall within sources. This is a result of the sta-
tistical independence of the source positions. Simple reasoning is sufficient
to show that the probability that a point chosen at random will lie within a
source is just the fraction of the total volume occupied by sources in the
vicinity of that point. Thus, we conclude that
W(r", Ar) = v2D(r") D(r" + Ar) , |Ar| > dimensions of v , (205)
and hence that
G(r',r") ~ m2D(r') D(r") , if |r' - r" | > dimensions of v . (206)
Hhen |r' - r" | has a value smaller than the source size, the determination
of W is slightly more difficult because the probabilities that r' and r" each
fall within sources are no longer independent. Instead, we have
(207)
where P....-(n,m) is the joint probability that U(r") has value n and that
simultaneously U(r" + Ar) has value m. By definition,
(208)
where Py,,(n|m) is the conditional probability that U(r") = n given that
U(r" + Ar) = m, and where Py.O) is the probability that U(r" + Ar) = 1.
From this definition and from Figure 33, we conclude after some reasoning
that
i , i ..i- ,,,, i \ n I r" \
PuU-Ollsr") = -v + (v "vvv } D(r } (209a)
where
V = (v - |Ar |)(v - |Ar |)(v - |Ar |) . (209b)
A A jr y i- t-
-------
117
To see this, note first that it is given that one point, say, r", lies in a
source, for example, the source marked v-, in Figure 33. Since the separa-
tion vector Ar is also given, the second point r" + Ar must lie somewhere
in the volume denoted v2 in Figure 33. If this second point happens to lie
in the shaded volume, i.e., v', then both U(r") and U(r" + Ar) will have
unit value; but if the point falls outside the shaded region, then U(r" + Ar)
will have zero value, unless the point happens to fall within a neighboring
source. The two events, r" + Ar lies in v' and r" + Ar lies outside v', are
mutually exclusive. Thus, Eq. (209) follows immediately.
FIGURE 33. RANGES OF THE VECTORS r" AND r" + Ar
GIVEN Ar AND GIVEN THAT r" LIES IN v
1
We found earlier that
= Pu(l;r" + Ar) = vD(r" + Ar)
but since |Ar| < v , v , v in the present analyses, we have
= vD(r")
(210)
Combining this result with Eqs. (209), (207), and (208), we obtain the informa-
tion required to evaluate W, given by Eq. (203), and subsequently G:
-------
118
G(r',r") =
-i
' - (v2 - vv1) D(r")J , |Ar| < dimensions of v .
(211)
where v' is given by Eq. (209b). In summary,
G(r'r") = -<
. * j
m D(r') D(r")
if any component of r1 - r" is larger than
the corresponding dimension of v ;
' -(v2 - vv1) D(r")] , otherwise;
(212)
where v1 = (YX -
v = vxvyvz
v - |Ar |)(vz - |Ar
Upon substituting this expression into Eq. (195), we can solve for the mean
~2
square inert concentration distribution Aj arising from the random array of
point sources.
For simplicity, let us assume that
D(r) = D = constant.
(213)
If the diffusivities K,, and K7 are also uniform in space, then all spatially
n ^, ~o ~
averaged quantities, such as A,, Aj, and S, will be independent of spatial varia-
bles. Consequently, regardless of the speed of the mean wind, there can be no net
advection of spatially averaged quantities, and we can therefore assume that u = 0
without any further loss of generality. In this case, we have
p(r-r',t,t') -
exp
H z
(x - x1)' + (y - yT (z - z')'
2 ?
2or," te
H z
(214a)
= 2KH(t - t1)
= 2K(t - t1)
(214b)
-------
119
The problem thus posed is one of a uniform spatial density of continuous
point sources acting in a homogeneous turbulence of arbitrary mean speed.
From Eq. (195), we now have
Aj(r,t) = m
•'/'(Iff
0 Ar >v
1/3
D p(r-r",t-t") p(r-r"-Ar,t-t') dAr dr" (215)
p(r_r",t-t")
Ar < v
1/3 v
p(r-r"-Ar,t-t') dAr dr"\ dt1 dt
where we have assumed that
= Vy = vz = v1/3
(216)
The Gaussian form [Eq. (214a)] of the kernels p permits a reduction of the
first integral within the braces in Eq. (215) to
1
1 -
^
v1/3
Ka^2 +
aS2)
1/2
L
1
A(c
v1/3 1
i1 + a" )
(217)
where
^ = 2KH(t - t1) , a[)2 = 2KH(t - t")
and so forth. The second integral within the braces in Eq. (215) reduces to
- 7vD)
(a'2 + a"2) (a'2 + a"2)1/2
(21;
-------
120
1 /3
In deriving Eqs. (217) and (218), we had to assume that a^, o^, » v . Fur-
thermore, we assumed earlier that vD «1. Consequently, Eq. (215) can be
written in the form
fa - *{(
D2 +
(2.)
3/2
2KH[(t -
t
D
')
+ (t -
t")]
3/2
(2K.
,)1/2
dt1 dt1
•2,A2 . Dm2t1/2
= m D t + o i /o—
- D
(219)
Similarly, we find that
72, ,x '2n2.2
Aj(r,t) = m D t
Making use of this result and Eq. (219), we obtain from Eq. (189)
(220)
,2/3
Dt3/2 1
t»
Y 2/3 „ 1/3
KH Kz
(221)
As indicated, Eq. (221) is not valid in the limit as t->0. This limiting
value is easily obtained, however, when one notes that at small times the
material is confined to a puff of volumes v equal to that of each source.
Straightforward analyses then lead to the result
(222)
The significance of these results to the behavior of A and B in the
problem under study becomes clear once Eq. (221) is incorporated into our
parameterization [Eq. (186)] and the result is used in Eq. (147). We obtain
SA •
-JJ-JT+ LA = mD - k'AB
(223)
-------
121
where
k- = kh(l + JW-T77) (224)
and
/2 -
a =
I/2
The parameter k1 is an effective rate "constant" whose departure from the
nominal value k is a measure of the SSV effects on the chemistry of A and
B. In the case where h - 1, we note that
lim k1 = k,
and that
k1 - Uj-^ . (225b)
The first of these indicates that SSV effects die out with time. From
Eq. (222), we find that
lim k1 - . (226)
In view of the second property, Eq. (225b) , the importance of SSV effects
diminishes as the source number density D, the turbulent diffusivity K, or
both increase.
Let us examine the SSV effects quantitatively. Because of the constancy
of Sn and SR, A is homogeneous in space. This fact results in
LA = 0
Assume also that SA = SB so that A = B. Equation (223) then reduces to
dt
(227)
-------
122
the solution of which is
A(t) =
mD
k1
- 1
(228)
where
Y =
It can be seen from Eqs. (228) and (225a) that the steady-state value of
A is independent of the SSV:
- limA(t) =
( 229)
Using A^, we can define a reaction time scale TR:
TR = (kAj-1 = (mkD)"'
(230)
This time scale is a rough measure of the time required for the chemical
reaction A+B->D to achieve a rate equal to that of the source emission rate.
Substituting TD into Eq. (224), we find that by the time the chemical decay
K
rate and the emission rate are in balance, the effective rate constant k1
has a value
a(mk)3/4 1
k1 = k
1 +
/4 K K 1/2
KHKz
(231)
Since k1 has a maximum at t = 0 and decreases thereafter, it follows from
Eq. (231) that if
mk
y = V4
D ' K,,
3/4
1/2
• 1
(232)
then SSV effects will play a significant role in the behavior of A and B as
they approach their steady-state values. This fact is illustrated graphically
in Figure 34, where we have plotted the ratio
A
(233)
-------
FIGURE 34. INFLUENCE OF SUBGRID-SCALE CONCENTRATION VARIATIONS
ON THE CHEMICAL RATE OF CHANGE OF .A
CO
-------
124
for several values of the dimensionless group y. Here, ASSV denotes the mean
value of A, which accounts explicitly for the SSV using k', and A denotes the
corresponding value obtained using k. According to Figure 34, y must have a
value comparable ti
are manifest in A.
q
value comparable to or larger than about 10 before significant SSV effects
It is of interest now to apply our criterion, Eq. (232), to an assess-
ment of the SSV effects in actual problems encountered in pollution simula-
tion studies. One such problem is the simulation of large urban areas
containing strong point sources, such as power plants and refineries. The
question here is whether such sources produce a significant SSV impact on
the predictions of an airshed model.
To approximate the conditions of a problem of this type, we let D
approach the value of one source per grid cell. In a typical case where
3 3
the grid cell dimensions are 10 x 10 x 50 m, we have
D = 2 x 10"8 m"3
We assume that the source emission rate in this case is 100 gm/sec. This is
a rough value for the NO emission rate of power plants. We also consider here
three conditions of atmospheric stability. These and the corresponding values
o
assumed for Ku and K are as follows (in m /sec):
n Z
Condition
Stable
Neutral
Unstable
KH
5
50
50
K
ID'1
1
25
Substituting the above values in Eq.' (232), we obtain the curves of y as a
function of k shown in Figure 35.
3
Note that y = 10 is the smallest value of y shown in the figure. Conse-
quently, all combinations of rate constant and stability classification shown
in this figure represent cases where significant SSV effects occur, at least in
-------
125
problems of the type considered here. For reference, we have indicated the
values of the rate constants for some of the important photochemical pollution
reactions. We should add, however, that for some of these reactions, the value
of y cannot be read directly from the figure because we assumed that the
concentration of the two reactants are approximately equal. In cases
where the reactant concentrations are in the ratio a (a
-------
«»•
^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
-------