Mathematical Modeling of Urban Air Pollution

General Theory
Robert G. Lamb and John H. Seinfeld’
Department of Chemical Engineering, California Institute of Technology, Pasadena, Calif. 91 109

The fundamental theory of urban air pollution mod- the pollutants resulted in alterations of the fluid behavior.
eling is developed. A number of points are considered in- From this point on we will not include the effect of tem-
cluding: The definition of the mean and random compo- perature on Ri.
nents of the wind velocity field, the necessary conditions For reasons t o be discussed later, it is necessary to rep-
for obtaining conventional Eulerian and Lagrangian mod- resent the wind velocity u = (u1, u2, us) in the form u =
els, the proper inclusion of nonlinear chemical reactions in il + u’, where it and u‘ are deterministic and random
a model, the nature of the spatial and temporal averages components, respectively. The deterministic component il
implicit in the computed concentrations, and the require- and the statistical properties of u‘ may be determined ei-
ments for simulating photochemical smog formation in an ther from observations or from a numerical simulation of
urban area. the fluid dynamics. We will subsequently consider a
means of determining these two components. Thus, Equa-
The development of urban air pollution models has tion 1.1becomes
been a subject of interest for a number of years. Unfortu-
nately, most of the published modeling studies have es-
chewed the theoretical problems associated with the prop-
er description of atmospheric diffusion and have concen-
trated instead on techniques for utilizing conventional
Gaussian plume and diffusion equation models. Before the Because urn’ is a random variable, Equation 1.2 repre-
capabilities of the various models and the meaning of sents a set of stochastic differential equations-i.e., corre-
their predictions can be assessed, four fundamental as- sponding to each sample function ua’ there is a e, satisfy-
pects of air pollution modeling must be considered: ing Equation 1.2. Thus, c, is itself a random variable,
The precise definition of the mean and random compo- which can be characterized by a probability density func-
nents of the wind velocity field that enters into a model tion. The probability density function for c L cannot be de-
The necessary conditions for obtaining various conven- termined in general, so we concentrate on the moments of
tional model forms from the general theory of turbulent the distribution of c,. The first moment of c L is its ensem-
diffusion ble average (c,(X,t)). Ensemble averaging of Equation 1.2
The proper way in which description of the rates of yields
chemical reactions among pollutants enters a model
The nature of the averaging implicit in the computed
concentrations-i.e., over what volume of space and period
W C J
of time are the concentrations averaged D , m + @,((el) + ~1’. ..., (e\) + et’)) + S,(X, t ) (L3j
We will consider each of these aspects, our objective
being to place within proper perspective the various ap- Let us consider the form of ( R , ( ( c l )+ CI’,. . ., ( c . ~ )+
proaches to air pollution modeling and to delineate those c q ’ ) ) . Suppose that among the N species M chemical
problems in need of further attention. reactions are occurring. For elementary reactions we may
The simulation of urban air pollution is essentially the express R , as
problem of describing the behavior of a number of chemi-
cally reactive species in the turbulent atmospheric bound-
ary layer. Consider N chemically reactive constituents in
a fluid. The concentration of each constituent must satis- where u,, is the stoichiometric coefficient for species i in
fy the continuity equation, reaction J, k , is the reaction rate constant for reaction j ,
and pi, is the reaction order of species 1 in reaction j .
In chemical reactions among air pollutants the presence
of very reactive species, such as free radicals, lead to some
reactions being very much faster than others. A common
R,(c,, ..., e,, T ) + S , ( X .t ) = 1, 2. ..., N (1.1) approach to problems of this type is the pseudosteady
where u,> is the tu-component of the wind velocity, D , is state approximation in which very reactive species, usual-
the molecular diffusivity of species i in air, R, is the rate ly present in low concentrations, are assumed to be
of formation or depletion of species i by chemical reaction, formed and destroyed a t approximately equal rates a t
T is the temperature ( a function of location and time), each instant of time. If it can be assumed that the induc-
and S, is the rate of injection of species i from sources. We tion period for the pseudosteady state is very short com-
follow the summation convention in which repeated Greek pared to the time scales of the molecular and turbulent
subscripts in a term indicate summation over the three transport terms in Equation 1.2, then for those species in
components of that term. Since pollutants occur in gener- pseudosteady state, R, = 0. This leads to algebraic equa-
al a t parts-per-million concentrations, their presence does tions for the pseudosteady state species in terms of the
not influence the meteorology. Consequently, the fluid ve- other species, which when ensemble-averaged yields
locities t~~~can be considered ipdependent of the c L . This
might not be true if a polluted layer were so concentrated
that absorption, reflection, and scattering of radiation by
where the number of species in pseudosteady state L < N .
ITo whom correspondence should be addressed If the total reaction order

Volume 7, Number 3, March 1973 253

 , such phenomenon as turbulent diffusion. Thus, turbulent
CP, 5
I
2
diffusion is an artifact of our lack of complete knowledge
then Equation 1.5 can be written as of the true velocity field. Consequently, one of the funda-
mental tasks in turbulent diffusion theory is to define the
deterministic and stochastic components of the velocity
field.
As we noted, u may be derived either from observations
z = 1,2, ..., L
or from the solution of a fluid dynamic simulation model.
These L equations are then coupled to the ( N - L )Equa- In either case the velocities are known only a t discrete
tions 1.3 for the species not in pseudosteady state (the points in space. This raises the question of whether these
accumulating species). In Equation 1.3, for the accumu- discrete values of u are sufficient to describe u at all
lating species points in space and time as required by Equation 1.1.
Assume that the three-dimensional velocity field is
(R((CI) + cl’, ..., (cn.) + CAV’)) =
sampled on a network of grid points with separations A X 1
= d l , AX2 = dg, and A X 3 = d3. The measured wind veloc-
ities can then be represented by u n m ( n l d l , ngdz, n3d3, t ) .
From these discrete measurements we wish to construct
Equations 1.3 and 1.6 now contain new dependent vari- u,(X1, Xz,X B , t ) . Let us assume that the measurements
ables, u a m are free of error. We can then invoke a three-dimen-
sional form of the sampling theorem (Papoulis, 1965;
Lamb, 1971b), a basic result in communications and sam-
so that now there are more unknowns than equations. pled-data control theory. According to this theorem, if the
This, in essence, is the closure problem associated with wind velocities do not have wave number components
the Eulerian description of turbulent diffusion. A number greater than n / d , , i = 1,2,3, then u,(X,t) can be com-
of approximate solutions to the closure problem have been pletely reconstructed from u a m ( n l d l , ngdz, n3&, t ) using
proposed (the principal ones are discussed in a paper by the relationship
Kraichnan, 1962) but each leads to an equation for the
concentration statistics which is accurate for only a limit-
ed class of problems. Among other approximate solutions (22)
is the gradient transport or mixing length hypothesis
which leads to the semiempirical diffusion equation (to be Let us explore the implications of Equation 2.2. To re-
discussed later). construct u(X,t ) exactly, the wave numbers associated
By contrast, in the Lagrangian approach, we attempt to with u must not exceed n / d l . Wave numbers are the in-
describe the concentration statistics in terms of the statis- verse of wavelengths, which can be thought of as the sizes
tical properties of the displacement of particles in the fluid. of turbulent eddies. In a typical urban area a lower limit
The difficulty in the Lagrangian approach lies not in on d , for weather stations is at best 1 km. Thus, in order
a closure problem but in determining the statistics of the for Equation 2.2 to hold, no turbulent eddies must exist
particle displacements. Also, this approach is not readily with sizes smaller than about 2 km ( 2 */Amln = x / d m l n ) .
applicable to situations involving nonlinear chemistry. However, we know that in the atmosphere there are
Nevertheless, because of the severity of the Eulerian clo- eddies as small as km. Instead of using uUmfrom
sure problem, we shall consider a Lagrangian formulation weather stations, we could perform a numerical simula-
as the basis for urban air pollution models. We begin with tion of the fluid dynamics. A typical mesh spacing for
a discussion of the determination of wind velocities. such a solution might be 0.1 km. In that case, Equation
2.2 would be valid as long as there are no eddies smaller
Definition of Deterministic and Stochastic than 0.2 km. Clearly, we would be moving in the right di-
Velocity Components rection as d r decreases, but if the smallest eddies are in-
The first problem with which we must deal is the defi- deed 10-4 km, we would need d , = 10-4 km. Such a fine
nition of the deterministic and stochastic velocity compo- grid, either for actual measurements or numerical solution
nents, &, and ua’. To illustrate the importance of this is obviously impossible. Thus, Equation 2.2 can only give
problem, suppose we have a puff of pollutant of known that portion of the velocity field with wave number less
concentration distribution c(X, to) a t time to. In the ab- than a i d , . That part of the velocity field containing wave
sence of chemical reaction and other sources, and if we numbers greater than a/d, must remain unknown and be
assume molecular diffusion to. be negligible, the concen- treated as random.
tration distribution c(X, t ) at some later time t is de- We therefore let La be the value computed from Equa-
tion 2.2. However, to use Equation 2.2 to determine a,,
scribed by the so-called advection equation:
we must filter out the components of u,m with wave num-
(2.1) bers greater than a i d , . This requires that we obtain infor-
mation on spatial scales of the turbulence a t a given time
If we solve Equation 2.1 with uu = k,, and compare the from continuous velocity records a t fixed locations. Unfor-
solution with observations we would find in reality that tunately, this cannot be done in general. Nevertheless, we
the material spreads out more than predicted. The extra need some simple way of estimating spatial scales of mo-
spreading is, in fact, what is referred to as turbulent dif- tion from time scales of motion in order to filter out the
fusion and results from the influence of the stochastic smaller spatial scales from unm.
component uu’which we have ignored. Now let us solve Perhaps the simplest approach to this problem is the
Equation 2.1 with the precise velocity field [Link] should following. First, we hypothesize that the time scale T[ of a
then find that the solution agrees exactly with the obser- perturbation in u L is roughly proportional to the spatial
vations (assuming, of course, that molecular diffusion is scale of the perturbation. Therefore, let us say that
negligible), implying that if we knew the velocity field
precisely at all locations and times, there would be no
254 Environmental Science & Technology
where K = ( K I , K Z , ~ 3 is) the wave number vector of the that the particle will not decay in the time interval (t’, t )
perturbation and is some characteristic velocity, say is

For such a particle, the probability of a transition from

where V is a volume with dimensions comparable to the X’ to X is simply the product of Q(X, t l X ’ , t’) and the
dI’s. The magnitude of the wave number vector, 1 ~ 1 , is probability that during the time interval it does not
roughly a measure of the inverse of the magnitude of the decay. Thus, in this case, Equation 3.1 becomes
spatial eddy sizes. Next, we examine the coherence func-
tions of the velocity records measured a t neighboring wind
stations. From these functions we estimate the period, T,,,,
‘3.2)
say, of the eddy with wave number components K~ = x / d I ,
i = 1, 2, or 3. Then, in view of the hypothesis (2.3), all we The density function, $, is relevant only when a single
need do is filter each velocity record (consisting of three particle is present. Suppose, however, that an arbitrary
velocity components) to remove all frequencies greater number n of particles are initially present and that the
than KIT,,,. position of the ith particle is given by the density function
Assuming Equation 2.3 is valid, once a time scale i I for $I(X, t ) . In this case it can be shown (Lamb, 1972b) that
each velocity component is determined, the data uam can the ensemble mean of the concentration c A v ( X , t ) aver-
be filtered by well-known techniques (Jenkins and Watts, aged over a sample volume AVcentered a t X is given by
1968) to produce a velocity record t&m, which satisfies the (3.31
conditions of Equation 2.2. Then, Equation 2.2 can be used
to give Q(X,t). Since we then know u(nld1, nzdz, n3d3, t ) , Hence, the so-called point concentration c(X, t ) , corre-
the true velocity, and Q(X,t), the deterministic velocity sponding to an infinitesimal sample volume AV, has the
(not just a t the measuring points but everywhere), the u‘ ensemble mean value
t,) = 2$,CX t )
are defined exactly a t the measuring points where they can
be analyzed for their statistical properties. Assuming the (a, i=1
(3.4)
statistics of the u‘ values are spatially smooth, we can in- By expressing the probability density qL(X’, t’) in
terpolate the statistics of U’ over the entire field. Equation 3.2 in terms of the known initial particle distri-
Having defined Q and u‘, we now return to the determi- bution (c(X0, t o ) ) and the known spatial-temporal distri-
nation of the mean concentration by Equation 1.3. In bution of particle sources S(X’, t’), say in units of parti-
spite of the fact that we now know both Q and the statis- cles per volume per time, and then substituting the re-
tics of u’, the determination of the statistics of cI from the sulting expression into Equation 3.4, we obtain the fol-
statistics of u’ is hindered by the closure problem. As we lowing general formula for the mean concentration:
have noted, there is a t present no generally valid solution (dX t ) )=
to this problem. We therefore turn to a Lagrangian ap-
proach in the next section where we assume that in place
of Q and the statistics of u‘, we have knowledge of Q and
the statistics of the Lagrangian velocity v’. Here v’ repre-
sents the velocity of a particle measured with respect to a
coordinate system moving with velocity Q starting from
This is the fundamental relationship for the mean con-
the point of release of the particle.
centration of a species which can decay exponentially in a
Equations Governing Mean Concentration of Particles in turbulent flow in which there are sources. The determina-
Turbulent F l o u tion of (c(X, t ) ) , given (c(X0, t o ) ) and S(X’, t ’ ) , rests on
the evaluation of the transition density Q(X, tlXo, t o ) . We
To avoid the closure problem associated with the Euler-
might note that if the turbulence is stationary and homo-
ian approach, we will derive the form of Lagrangian mod-
els for turbulent diffusion and chemical reaction. We geneous Q(X, t l X o , t o ) = Q(X - XO; t - to)-i.e., the
transition density depends only on the displacements in
begin by considering a single particle which is a t location
time and space and not on where the particle is or when it
X’ a t time t’. The subsequent motion of the particle can
was introduced into the flow. In the case, empirical data
be described by its trajectory, X[X’, t’; t ] . Let $ ( X , t ) be
indicate that the transition density Q(X - XO; t - t o ) is
the probability density function for finding the particle
marginally and perhaps jointly normal (Batchelor, 1949),
within an infinitesimal volume centered a t X a t time t.
For generality, we can assume that the initial position of
the particle is not known precisely but rather is given in
the form of a probability density function, $(X’, t’). The where CT is the transpose of the column vector ( which
probability of finding the particle a t X a t t is simply the has elements
product of the probability that the particle undergoes a
displacement from X’ to X in (t’, t ) and the probability c = x - x , - ( X - X , ) 1 = 1 2,3 ( 37 )
that the particle was a t X’ a t t’, integrated over all X’. If and where P-I and IPJ are the inverse and determinant,
Q(X,tl X’,t’) is the conditional probability that if the par- respectively, of the covariance tensor P,
ticle was a t X‘ a t t’ it is a t X a t t, then p / = KL,) (3 8)

$CX. t ) = J‘J=J=Q(X, tl X . t’)$(X’. t’)dX’ ( 3 1) If the coordinate axes are the principal axes of the tensor
-a -_ -r
P, then PI, = 0 for 1 # J We assume that the mean dis-
If we assume, in addition, that the particle may decay placements ( X I - X , O ) are due only to the deterministic
according to a first-order rate law, the “lifetime” of an in- velocity
dividual particle is exponentially distributed. If the decay (X,- x,) = J [ i i ( X , t’Mt’ ( 39)
coefficient is a function of time, A ( t ) ,then the probability

Volume 7, Number 3 , March 1973 255

where 8’ is the particle position at t’ if the only velocity possible to determine either experimentally or theoretical-
were iil. ly in nonstationary, inhomogeneous turbulence. Perhaps
If P,, = 0 for i # j , Equation 3.6 becomes, for T = t - the only feasible way of determining Q is by numerical
to, simulation of turbulent flows, such as reported by Lamb
QCX - xo;t - t o )= (1971a) who used the fluid simulation model of Deardorff
(1970).
When nonlinear chemical reactions are occurring the
Gaussian puff and plume formulas cannot be used, al-
If, in addition, there is no mean flow, ( X I - X , O )= 0 , and though for very slow reactions Friedlander and Seinfeld
rl is replaced by ( X , - X , O ) .We often denote the diagonal (1969) have extended the Lagrangian similarity theory of
elements of P, P I l ( 7 ) ,by u 1 2 ( ~ )If. the turbulence is also Batchelor to enable inclusion of such reactions. In general,
isotropic, there is no preferred direction and C T , ~ ( T )= however, these formulas cannot be used as the basis of a
fJ2(T). general urban air pollution model. Let us return to Equation
Let us consider some special cases of Equations 3.5 and 3.5 in a n effort to develop a general Lagrangian model.
3.6: One of the disadvantages of using Equation 3.5 in prac-
(a) Instantaneous point source at XO which emits unit tice is that it requires a tremendous amount of computa-
concentration of pollutant with a reflecting barrier at X s tional effort if the time interval, t - t o , is large and there
= 0 (the ground). No decay of pollutant. are many particle sources. This is explained by the fact
MX, t)) = that Equation 3.5 essentially treats each particle individ-
ually from the time that the particle is produced to the
time that the calculations end. The computational effort
is rewarded, however, by the achievement of high resolu-
tion descriptions of the spatial and temporal variations in
the particle concentration statistics. If one is willing to
sacrifice part of this resolution or if a fluid simulation
If, in addition, there is a mean wind 21 = U in the direc- model is used to reduce the scales of the turbulent veloci-
tion of the X I axis, the first argument of the first expo- ty components, the computing time required by Equation
nential in Equation 3.11 becomes - ( X I - X I O - UT)’/2 3.5 can be reduced. Let us see how this might be done.
u12(7). Consider the situation in which the random component
(b) Continuous point source of unit strength a t XOwith of the velocity of any particle, u , ( t ) . has a correlation
iL1 = U, 22 = iL3 = 0. Reflecting barrier at X 3 = 0. No decay function R,,(t; T ) = ( u l ( t ) u , ( t + 7 ) ) which vanishes suffi-
of pollutant. ciently rapidly with increasing 7 that a time scale
n-

exists for all values of t and all possible points (X’, t’) of
release of the particle. Under this condition it is to be ex-
pected that the motion of any particle a t any time t will
be statistically independent of its motions prior to the
time t - At where
-
Since the source is continuous we want the solution as t
a. At steady state, there is no time dependence so we
must convert from travel time 7 to distance X I - X l o .
At >> maxT,,
J. J
(3.15)

-
Then, the u12 becomes a function of distance from the
source, X I - X l o . Letting t m in Equation 3.12 gives
In particular, given a sequence of times t l , tz, . . . t , =
n A t , where At satisfies Equation 3.15, we may expect that
the conditional density function Q(Xn, tnlXn-1, t n - 1 , . . . ,
XO, t o ) which gives the probability t h a t the particle is a t
Xn at t , given that it was at Xn-1 a t tn-1, and so forth,
satisfies

Q(X,,,t,,IXn-i. t n - ) . ..., Xo. to) = QcX,,. tnIX,,-i,t,,-ij (3.16)

This result is the so-called Gaussian plume equation for a
source of unit strength. In this case the particle displacements constitute a Mar-
The formulas in (a) and ( b ) are the classic relationships kov process and Q(X,, tnlXs, t s ) satisfies the Chapman-
for diffusion in homogeneous turbulence. The continuous Kolmogorov equation
Gaussian plume formulas have been used extensively in
predicting pollutant dispersion from point and line
sources. They have been fairly successful in predicting
Q(X,,.[Link],. t.) = Jc
-,
Jr
- ?
s=Q(X,. t,lX,, t,“@(X,. t,,[Link],
-1

(3.17)
long-time averages of concentrations near the source (XI n>m>s
- XI0 < 10 km) under steady meteorological conditions
when used in conjunction with empirically determined Noting that Q(X, t l X ’ , t ’ ) in Equation 3.5 can be replaced
values of the u12, such as from the Pasquill-Gifford curves by Q(X,, tnlX,, t s ) , we can combine Equations 3.5 and
(Pasquill, 1962; Slade, 1968; Turner, 1969). Although rig- 3.17 to yield
orously these formulas do not apply in nonstationary,
inhomogeneous turbulence, they can provide reasonable l. t,,-,)) x
MX,, t,,)) = J : S _ : S _ a ~ ( ~ , . t , ~ ~ ; , -t,,-i)(e(xr,-l.
order of magnitude estimates in many practical circum-
stances. These formulas are used not only for their sim-
plicity but also because Q(X, tlXo, t o ) is practically im-
256 Environmental Science & Technology
 exp( - J r A ( ~ ) d ~ ) d t ~ d X(3.18)
,
exp( - A(t)At) = 1 - A(t)ilt f o(b*) (3.30)
Equation 3.18 is a general relation for the mean concen-
tration of exponentially decaying particles under the re-
striction that A t is significantly greater than the Lagran-
gian time scale of the turbulence. We would like to reduce
Equation 3.18 to the differential equation forms for (c(X,
t ) ) more commonly encountered. T o do this, we must Substituting Equations 3.27-3.31 into 3.26, we obtain
make additional assumptions. First, it is necessary to as-
+ O(At*)= -(c)A(t)At a +
- [l - A(t)~lt]d~,,((e)(AX,,))
sume that temporal variations in S(X, t ) and A ( t ) are at
small compared to At and that spatial variations in S(X,
t),are gradual, namely ;(i - A(t)At)-
2 dX"dXd
+
((cXAX,,AXd)) S(X tjAt -
(3.19) a (S(AX,,))At+ O((AXAX,jAX-,)) + Or(Ax,AX,)At) (332)
dx,,
We must now evaluate (AX,) and (AXdX8). The mean
displacement in a time interval A t is
(AX,,) = iC,,(X, t)At (3.33)
Secondly, we must assume that the spatial and temporal The mean square displacement tensor is defined by
inhomogeneities in the turbulence are of such scales that
the transition density, Q, of a particle released anywhere in (AX,,AXj ) = smJEs'AX.2AX jQ(AXiX)d(AX) (3.34)
- 7 -_ -r

the fluid is the Gaussian density for travel times t in the

On the basis of the assumption that Q is Gaussian and
range 0 I t I A t . Actually, this condition is probably au-
using Equation 3.33, we can write
tomatically satisfied when Equation 3.15 is true. Thus,
over each interval of A t a particle behaves as though it (AX,,AX,j) = K,,j(X, t)At + k ( X , f)G,j(X,tMt' (3.35)
were in a field of stationary homogeneous turbulence. We
can refer to such a situation as locally stationary and ho- where K<,J= ( L ' , , ~ J ) T+, ,(, ~L ~ C , ) T ~ ~ ~ (3.36)
mogeneous turbulence. Quantitatively, the fluid velocity Substituting Equations 3.33 and 3.35 into 3.32, dividing
components should satisfy the following conditions: by A t , and letting At 0, we obtain
+

(3.23) which is recognized as the Fokker-Planck equation corre-

sponding to the random process represented by Equation
3.26 (Chandrasekhar, 1954). We should point out that to
achieve the limit At 0, we must consider At (and, for that
+

matter, all time variables in Equation 3.32) as having been

nondimensionalized by an arbitrary time scale T. The limit
A t 0 is then achieved by letting T + m .
+

If the coordinate axes coincide with t h e principal axes

Under these conditions the transition density Q(X,, of the Ka3 tensor, or if the turbulence is isotropic, then Kea
= 0, LY # p, and Equation 3.37 reduces to the customary
t n l X n - l , t n - l ) can be written Q(AXjX), the transition
probability that the particle will be displaced a distance diffusion equation
AX from X during At. In addition, this transition density
will be Gaussian.
Under all of the above conditions, Equation 3.18 takes
the form When we go one step further and assume that K,, is

(dX. t + At)) = s JzJr

- -_ -
@ANX - AX) exp( - A ( t ) A t )X
spherically symmetric and constant throughout space, we
find that Equation 3.38 reduces to

[(dX - AX. t ) ) + AtS(X - AX, tj]d(AX) (3.26)

The objective is to derive a differential equation for (c(X, where

t ) )from Equation 3.26.
Expanding (c(X, t + A t ) ) , Q(AXjX - AX), (c(X - AX, In spite of the fact that Equations 3.37-3.39 are ex-
t ) ) , exp ( - A ( t ) A t )and S ( X - AX, t ) in terms of Taylor pressed in an Eulerian coordinate frame, they are based
series about the point (X, t ) ,we obtain on Lagrangian statistics, as can be seen from the defini-
tion of the "diffusivity" Kn3. From a strictly mathemati-
(dX.t + At)) + 7ace)
i A t + O(At2)
(dX. t ) ) (3.27) cal point of view, Equations 3.37-3.39 are valid only for
infinitely large time scales T. From a practical standpoint,
aQ
ax AXck+
Q(AXIX - AX) = Q(AXIX) - __
t
however, it is necessary only that T >> A t . Also, we have al-
ready seen that A t >> maxT,, so we conclude that
!ax,,ax, AX,,AXi - ... (3.288) maxT,, << At << T
,
(3.40)

Volume 7, Number 3 , March 1973 257

 The basic Lagrangian Equation 3.5 applies only to linearly An eddy diffusivity tensor K,B exists, defined by
reactive particle species, because in deriving the funda-
mental relationship, Equation 3.2, on which this equation (4.1)
is based, we assumed that the chemical reactions occur and that if reactions are occurring, moments of the form
independently of the particle displacements. This as- (ci’cj’) can be neglected. Generally it is assumed that it,, =
sumption is clearly not always valid for nonlinear reac- 0 for CY + p. Then Equation 1.3 becomes
tions where the reaction speed is determined by the fre-
quency of particle collisions. To handle this case using
Lagrangian statistics is very complicated (Lamb, 1972c),
but it appears that this approach may be the only way a t
present whereby the complex and important microscale
transport and reaction phenomena which occur in the vi-
It can be shown (Lamb, 1972a) that the two above as-
cinity of point and line sources can be described. How-
sumptions, and hence Equation 4.2, are probably valid
ever, the Lagrangian theory for nonlinear processes ap-
provided that in addition to Equations 3.19, 3.21, and
pears to be too complex to provide a viable basis for mod-
3.22-3.29, the following condition is met: that the time
eling macroscale features of the concentration distribution
scale of the fastest reaction described by R(c1, . . . C N ) be
in urban atmospheres. For this job any one of the Equa-
much larger than the Lagrangian time scale, T , of the tur-
tions 3.37-3.39 may be used. bulence. This list of restrictions is essentially identical to
Regarding these equations,’we point out that each may be
that which was placed on Equation 3.38, so as far as
applied to certain nonlinearly reactive substances owing to
applicability is concerned, Equations 3.38 and 4.2 are
the assumptions in Equations 3.19-3.21 on which each of
identical. The only significant difference between these
these equations rests.
two equations lies in their diffusion terms. In Equation
By virtue of these assumptions the temporal and
3.38, K,, is dependent strictly on the Lagrangian statistics
spatial variations in the mean concentration field are so
of the turbulence while in Equation 4.2, K,, is an Euler-
much larger in scale than the Lagrangian time scale T =
ian parameter-Le., a function of each point in space,
max Tl, and the average distance L = [(a,z + ( u 1 2 ) ) ~ z ] 1 / 2
which is to be determined by empirical data. However,
that a particle travels in time T that the collision frequen-
since the Lagrangian properties on which the K,, are de-
cy, and hence nonlinear reaction rate, experienced by any
pendent are very difficult to measure directly, the K,, are
particle during any period T , is independent of the parti-
usually determined by using them as parameters for fit-
cle’s path and moreover is just that frequency, or reaction
ting the solution of Equation 3.38 to empirical concentra-
rate, which is appropriate for the local mean concentra-
tion data. This same procedure is used to determine the
tion. Thus, we may use the terms A(c,) and S,to repre-
sent not only the linear reaction rate and strength of the
K,, so that all in all !here is little significant difference
between Equations 3.38 and 4.2. Because of the empirical
systematic sources, respectively, of the ith species, but
nature of the diffusivities, the accuracy of both equations
also the sources and sinks resulting from the nonlinear
is also dependent on the degree to which the conditions of
chemical reactions among the N species. Specifically, the
the problem to which these equations are applied corre-
term A(c,) can represent all those reactions which deplete
spond to the conditions under which the diffusivities were
species i, . ..l .. measured.
As we have already pointed out, Equations 4.2 and
3.37-3.39 apply only to time and space scales which are
and similarly S,(X, t) can represent the formation of much larger than the corresponding scales T and L of the
species i by reactions among the other species present, turbulence. Thus, to regard these as so-called point equa-
M“ tions is to imply that a “point” has spatial dimensions
S,(X t ) = ~VLJk,i?i(Cl)”J. (3.42) much larger than Ll(i = 1, 2, 3) and that an “instant” is a
,=I 1-1
ICt period long compared to T . Let these minimum resolvable
It turns out that the values of 7 and L, characteristic of scales be denoted by 1, and T-Le., I , >> L,, T >> T . Now in
the atmosphere, are so large that conditions in Equations applications where the chemistry is nonlinear, analytic so-
3.19-3.21 are grossly violated near isolated sources such as lutions are not available for these equations so they must
smoke stacks and highways. As a result Equations 3.37- in general be solved numerically. Depending on the type
3.39 provide only an approximate description of the ma- of computer available, it may prove infeasible to select a
croscale features of the mean concentration distribution difference grid with space increments AX, = 1, and a time
because they completely ignore the effects of the micros- increment A t = T. When intervals AX, > 1, and/or At > T
cale phenomena (mentioned earlier) induced by such are necessary, an additional problem arises in connection
sources. Since these microscale phenomena can have a with the nonlinear terms.
profound effect upon the large-scale concentration fea- To illustrate the nature of this problem, we consider for
tures in cases where nonlinear reactions are active, we convenience the one-dimensional equation
may expect that Equations 3.37-3.39 will prove to be in-
adequate for simulating chemically reacting air pollution
until some method is devised for parameterizing the ef- with K x a constant.
fects of the microscale phenomena in terms of the large Suppose we solve this equation numerically using At =
scale concentration distribution. This is one of the impor- T and AX > lx. Since a discrete grid with spacing AX
tant problems of air pollution modeling which remains to cannot resolve features of ( c ( X , t)) which have wave num-
be solved. ber components K > a/AX (by virtue of the sampling the-
orem mentioned earlier), the difference equation represen-
Semiempirical Equation of Atmospheric Diffusion tation of Equation 4.3 which we solve numerically can de-
The most common approach to problems of atmospheric scribe only the variations in (c(X, t ) ) which have wave
diffusion has been based on the Eulerian description of number scales K I n / A X . Moreover, the difference equa-
Equation 1.3, together with two assumptions that tion cannot account for the effects of perturbation in ei-

258 Environmental Science & Technology

ther (c) or S with K > a / A X . One way of remedying this Simulation of Chemically Reacting Air Pollution
situation is to filter from both (c) and S all wave numbers
greater than a / A X . This can be done by space averaging In the previous sections we have developed both La-
as follows: grangian and Eulerian descriptions of air pollutant behav-
z 1 X+AX
ior.
( c ( X ,t ) ) = (C(X’, t))dX’ (4.4) The question now arises-what do the restrictions on
2AXS,_,,
~

the use of the Lagrangian and Eulerian models mean in

with a similar definition for S ( X , t ) . If a(K, t ) and b ( ~t ,) practice? T o answer this question, we consider the prob-
are the Fourier transforms of (c) and S, respectively-i.e., lem of simulating photochemical smog, which is undoubt-
edly the most challenging current problem in air pollution
( d X ,t ) ) = ~ u ( Kt ) ,e i K y dand
~ . so forth ( 4.51
modeling. The first step is the formulation of a kinetic
then it is easy to show that the corresponding transforms mechanism for photochemical smog-Le., the determina-
of ( E ) and S are C ( K , t ) = (sin KAX/KAX)U(K, t ) , b ( ~ t,) = tion of a functional form for R1.
(sin K A X / K A X ) ~ t( )K. Due
, to the nature of the weighting The nature and characteristics of atmospheric contami-
function sin KAx/KAx, fi and 6 are virtually zero for nants suggest certain difficulties in the formulation of a
values of K > */AX, and thus ( E ) and S have properties kinetic mechanism of general validity. First, there is mul-
which are compatible with our difference equation version tiplicity of stable chemical species in the atmosphere.
of Equation 4.3. Our aim, then, is to use the difference Most species are present a t very low concentrations,
equation to obtain the space averaged concentration (E). thereby creating major problems in detection and analy-
T o obtain the differential equation for ( E ) we could av- sis. A number of atmospheric constituents, in fact, proba-
erage Equation 4.3 in the manner of Equation 4.4, but bly remain unidentified. Second, there are a large number
this task is more easily performed in wave number space. of short-lived intermediate species and free radicals which
We substitute Equation 4.5 and the corresponding rela- participate in an enormous number of individual chemical
tionship between S and b ( ~ t,) into Equation 4.3 and ob- reactions. However, while we must admit t o only a partial
tain for the case of k = 0, understanding of atmospheric reaction processes, it re-
mains essential that we attempt to formulate quantitative
descriptions of these processes suitable for inclusion in an
overall simulation model.
Multiplying by sin K L X I K A Xwe get A suitable mechanism must not be overly complex,
since computation times for the integration of the basic
model, in which the mechanism is to be imbedded, are
likely to be excessive. On the other hand, too simplified a
mechanism may cmit important reaction features. A
the inverse of which yields major issue in this regard is that the mechanisms predict
the behavior of a complex mixture of many hvdrocarbons,
(4.7) yet do so with a paucity of detail. The goal, then, is to
achieve acceptable accuracy in prediction without an
In other words, for linear problems ( E ) satisfies an equa- undue computational burden.
tion identical to that governing ( c ) with the forcing term Kinetic mechanisms that have been proposed fall into
in the former case being the space averaged source distri- two general categories: detailed mechanisms for photoox-
bution. idation of a single hydrocarbon and compact generalized
The situation is entirely different, however, when nonlin- mechanisms for a complex mixture. Detailed mechanisms
ear terms are present in Equation 4.3. Placing Equation which attempt to account for the history of all species
4.5 and the analogous expression for S into Equation 4.3, generated must be ruled out for three reasons. First, while
we get the aim of those developed thus far has been completeness
of description, this thoroughness has been achieved
au(K , t -
~-
at - K , K ~ c I (tK) ,- k
s U(K - K’)u(K’)~K’
It is immediately clear from this equation that in contrast
+ b ( ~t ), (4.8) through the inclusion of a number of reaction steps that
involve free radicals. Unfortunately, knowledge of the
rates of these reactions is imprecise. Furthermore, when
to the linear problem, the low wave number components several free radical reactions are included in a mecha-
of ( c ) which comprise ( E ) are affected, through the nonli- nism, the flexibility in the choice of rate constants is in-
near term, by all wave number components of ( e ) . This creased, as each imprecisely known parameter can be var-
presents the difficult problem of representing the effect on ied independently in the process of matching prediction
( E ) of the high wave number components of ( c ) which lie and experiment. To the extent that detailed mechanisms
outside the spectrum of ( E ) and are unresolvable by the possess this flexibility in parameterization, the validity of
difference grid. There is no rigorous solution to this prob- comparison of prediction and experiment is diminished.
lem because there is no exact way of expressing the effects Second, computation time is a limiting factor in the solu-
on ( E ) of subgrid scale variations in S in terms of and s tion of the coupled partial differential equations that
(2.). Perhaps by solving the concentration equation in wave comprise the overall airshed model. The inclusion of a
number space (Equation 4.8 in the present instance), we detailed mechanism in such a model greatly increases the
could treat a larger range of the spectral components of the computational burden and is to be avoided if a t all possi-
concentration and source distribution explicitly and thereby ble.
reduce to a negligible level the effects of the subgrid scale Finally, the decision to develop and implement a de-
or truncated components. Obviously, one can conceive of tailed mechanism implies the desire to represent reaction
problems in which it would be infeasible to eliminate the processes as accurately as is feasible. Thus, a relatively
subgrid scale effects by this method. What is needed in the large number of reaction steps must be incorporated in the
long run is a suitable heuristic parameterization of the ef- description of the dynamics of consumption of a particular
fects of the unresolvable fluctuations. The development of hydrocarbon, such as propylene. Reaction dynamics will,
such a parameterization will be the topic of a future study. however, vary for the many hydrocarbon species present
Volume 7, Number 3, March 1973 259
 ~~ ~~~~~

in the atmosphere. If, for example, 30-40 steps are re-

quired to describe propylene kinetics, and 50 hydrocarbon Table I . Generalized Mechanism for Photochemical Smog
species, each having unique dynamics, are believed to Reactions Rate constants, ppm, min
exert a significant impact on atmospheric reaction pro- 1. N O z + h u + N O + O 0.355 min-
cesses, one is faced with an intractable representation of 2. O + O 2 + M - O 3 + M 2.76 X 1 0 6 m i n - l a
the system. 3. 0 3 +
+
NO
-
NOz O2
-+ +
+
21.8 ppm-’ min-’

1
To summarize, then, we require a mechanism which: 4a. 0 3 NO2 NO3 0 2
describes reaction rate phenomena accurately over a spec- 4b. No3 NO2 N205 6 X p p m - l min-’
ified range of concentrations, is a parsimonious represen- +
4c. N205 H20-2HN03
tation of the actual atmospheric chemistry, in the interest + +
5. NO NO2 H 2 0 - 2 H N O z 2.5 X 1 0 - 3 p p m - 1 m i n - l b
of minimizing computation time, and can be written for +
6. HNO2 h~ -L NO OH. + 5 x 10-3min-’
general hydrocarbon species, with the inclusion of variable
stoichiometric coefficients, to permit simulation of the
behavior of the complex hydrocarbon mixture that actual-
7. CO
8. HOz.
+ OH- 4 COz
+
NO2
9. H C + O + a R O 2 *
-
HNOz
+
HOy
0 2 +
2 x 1O’ppm-l min-’
1 X 10’ p p m - ’ min-’
3.1 X l o 4 ppm-’ min
( a = 5)
ly exists in the atmosphere. In short, we seek a mecha- 10. H C + O ~ - P R O ~ . + Y R C H O 1 . 7 X 1 0 - 2 p p m - 1 min-’
nism which incorporates a balance between accuracy of ( p = 1.9)
prediction and ease of computation. 11. HC + O H v - 6 R 0 2 . + t RCHO 1 X 1 0 4 p p m - ’ min-’
The two mechanisms which fulfill the above require-
ments for suitability and have been most successful in
simulating smog chamber data are those of Eschenroeder
12. ROz* + NO --
NO2 + 0 OH*
+ NO2 PRODUCTS
13. R 0 2 -
(6 = 0.2; t = 0.22)
1.8 X l o 3 ppm-’ min-’
10. ppm-’ min-’
and Martinez (1972) and Hecht and Seinfeld (1972). The 14. HO2. + NO NO2 + OH.
--+ 1.8 X 1 0 3 p p m - ’ rnin-’
Hecht-Seinfeld mechanism is reproduced in Table I.
NOTE: It is possible to combine Reactions 4a-4c as
We might note a few points concerning this mechanism.
All reactive hydrocarbons have been lumped into one fic- 03+NOz 2HN03
H20
titious species, called HC. Also, all of the peroxy radicals
in which the overall rate is that of Reaction 4a, the rate-controlling step.
(those capable of oxidizing NO to NO2) are lumped into Pseudofirst order.
the species Ron.. Note that R 0 2 . is merely symbolic of Pseudosecond order.
these radicals; some may not contain exactly two oxygen
atoms. In reactions 9-11, the stoichiometric coefficients a,
(3, and 6 govern the growth of the total peroxy radical pop-
the i, j element of which isdR,/dc, and are the n solu-
ulation, whereas 9 and e describe the amount of aldehyde
tions of the equation, det(J - M) = 0.1
production.
First we consider a measure of the time scale for the
If CO and H2O are present, we treat HO2. and ROz.
fastest reactions in the system. Perhaps the most charac-
as separate species; in the absence of CO and HzO, both
teristic set of fast reactions in the smog system are reac-
are included in Ron.. In the latter case the coefficient, 8,
tions 1-3, the principal inorganic reactions. Consider the
cannot be interpreted exactly as the fraction of RO2.
dynamics of a system initially consisting of pure NO2 and
which is HOz.. This is because peroxy radicals other than
in which only reactions 1-3 are occurring. The system is
HOz. can, ultimately be degraded to OH.. For example, described by three differential equations (one each for
the history of a typical acyl (RCO) radical arising from NO2, 0, and Os). At t = 0 the largest nonzero eigenvalue
atomic oxygen attack on an olefin might be is -2 X 106 min-1, the pseudofirst-order rate constant for
Reaction 2 in Table I. Thus, the time required for this
0
11 0
RC. 4 RCOO.
0
11 -NO t o NO,
0
1I
RCD R. 0
system to reach steady state is the order of 10-6 min.
Clearly, the steady state approximation for oxygen atoms
is valid. Similar calculations will confirm the steady state
approximation for the other radical species in Table I.
Now consider the full mechanism in Table I. Our objec-
tive is to determine under what conditions the reaction
rates, R N O ~RNO,, R o ~and RHC, corresponding to this
mechanism and the parameters in Table I can be included
in Equation 3.38, the basic Lagrangian model. To include
Thus, an original acyl radical participates in the conver- these reaction rates, perhaps split into -A(c,) S((cl), +
sion of three molecules of NO to NO2 with an OH radical . . . , (cN)),the time scale for the system must satisfy Equa-
resulting. tions 3.19 and 3.20. As we noted, a measure of the time
The mechanism in Table I has four accumulating scale of the system is its eigenvalues. For the mechanism
species: NO, NO2, HC, and 0 3 . The following species are in Table I, the four eigenvalues a t three values of time are
in pseudosteady stage: 0, OH., HOz., ROz., “ 0 2 , given in Table 11. The system is characterized by a wide
H N 0 3 , and RCHO are products, while CO and H20 are separation in the magnitudes of the eigenvalues. The large
generally present in great excess and can be considered of negative eigenvalue corresponds to processes taking place
constant concentration. The dynamics of a batch system on time scales much faster than the rate determining
will then be described by four differential equations. steps. (It is important to note that in spite of the fact that
Earlier, we developed relative criteria for validity of the we have employed the steady state approximation for free
Lagrangian models based on time scales. Our main con- radicals, we still have a wide separation in the A’s.) After
cern with the photochemical smog mechanism is the time 1 min, the component of the linearized solutions corre-
scales associated with the system. A measure of the time sponding to X I will have a magnitude the order of e-30,
scale for the chemical reaction system is given by the and will be negligible. Thus, the important rate-determin-
eigenvalues of the reaction rate equations. [The eigenval- ing eigenvalue is X2, which is of the order of 10-2.
ues Xi, i = 1, 2, . . . , n, of a system of n ordinary differen- In considering Equation 3.38, A2 assumes the role of
tial equations are those of the n X n Jacobian matrix J , A ( t ) . The condition in Equation 3.20 implies, then, that

260 Environmental Science & Technology

Table II. Eigenvalues of Mechanism in Table I

t, rnin A1 A2 A3 A4

1.3 -33.1 4.63 X 1.24 x 10-4 1.28 x 10-7

60.3 -21.2 3.51 X 4.83 x 10-4 - 6.59 x 10-4i 4.83 x 10-4 + 6.59 x 10-4i
119 -21.5 -8.98 x 10-3 -2.04 X l o - ’ - 1.71 X 10-2i -2.04 X + 1.71 X 10-2i
NOTE: Initial conditions used were CNO,(O) = 0.04,CNO(O) = 1.5. c 0 3 ( 0 ) = 0,CHC(O) = 3.0.

At must be the order of 10 min or less. Fortunately, this Jenkins, G. M., Watts, D. G., “Spectral Analysis and Its Appli-
estimate is consistent with Equation 3.40, since max Tij is cations,” Holden-Day, San Francisco, Calif., 1968.
Kraichnan, R. H., Symposia i n Applied Math., 13,199 (1962).
generally the order of 1 min. This estimate of A ( t ) corre- Lamb, R. G., “Application of a Generalized Lagrangian Diffusion
sponds to the initial conditions in Table 11. These were Model to Urban Air Pollution Simulation Studies,” Proceedings
chosen as typical of those in a smoggy atmosphere. of the Symposium on Air Pollution, Turbulence and Diffusion,
Therefore, we can conclude that Equation 3.38 will rep- Las Cruces, N.M., 1971a.
resent a valid urban air pollution model for photochemical Lamb, R. G., “Note on the Application of K Theory to Diffusion
Problems Involving Nonlinear Chemical Reaction,” A tmos. E n -
smog as long as the concentration field is smooth (Equa- uiron., in press (1972a).
tion 3.21) and the concentrations are a t levels such that Lamb, R. G., “Numerical Modeling of Urban Air Pollution,” PhD
the rate-determining eigenvalue allows us to use a A t of thesis, University of [Link] t Los Angeles, Calif., 1971b.
the order of 10 min. Lamb, R. G., “Statistics of Marked Particle Concentrations in
Turbulent Fluid, Part 1: General Theory,” J. Fluid Mech., sub-
T o meet these requirements it will be necessary to mitted for publication, 1972b.
space-average the source inputs to remove steep gradients Lamb, R. G., “Statistics of Marked Particle Concentrations in
in concentration as would result, for example, from free- Turbulent Fluid, Part 2: The Case of Chemically Reactive Par-
ways and power plants. This point has been mentioned. We ticles.” ibid., submitted for publication, 1972c.
reiterate that the parameterization of strong localized Papoulis, A., “Probability, Random Variables and Stochastic
Processes,” McGraw-Hill, Ncw York, N.Y., 1965.
sources of reactive material for inclusion in a large simula- Pasquill, F., “Atmospheric Diffusion,” D. van Nostrand, London,
tion is still a n important unsolved problem. 1962.
A description of the simulation of photochemical smog Roth, P. M., Reynolds, S. D., Roberts, P. J . W., Seinfeld, J. H.,
in the Los Angeles basin by means of numerical solution “Development of a Simulation Model for Estimating Ground
Level Concentrations of Photochemical Pollutants,” Systems
of the semiempirical equations of atmospheric diffusion is Applications Inc., Beverly Hills, Calif., 1971.
given by Roth et al. (1971) and Seinfeld et al. (1972). Seinfeld, J. H., Reynolds, S. D., Roth, P. M., “Simulation of
Urban Air Pollution,”Adu. in Chem. Ser., No. 113 (1973).
Literature Cited Slade, D. H., “Meteorology arid Atomic Energy 1968,” [Link].
Energy Comm., TID-24190, J u l y 1968.
Batchelor, G. K., Aust. J. Sci. Res., 2, 437 (1949). Turner, D. B., “Workbook of Atmospheric Dispersion Estimates,”
Chandrasekhar, S., “Stochastic Problems in Physics and Astron- U.S. Dept. of Health, Education and Welfare, 1969.
omy,” in Selected Papers on Noise and Stochastic Processes, N.
Wax, Ed., Dover, New York, N. Y., 1954. Receiued for review May 3, 1972.. Accepted October 4, 1972. Pre-
Deardorff, J. W., Geophys. Fluid Dynam., 1,377 (1970). sented at the Symposium on Science in the Control of Smog,
Eschenroeder, A. Q., Martinez, J. R., “Concepts and Applica- California Institute of Technology, Pasadena, Calif., November
tions of Photochemical Smog Models,” Adu. i n Chem. Ser., No. 1971. Work performed in part at the University of California, Los
113 (1973). Angeles, with Federal funds from the Environmental Protection
Friedlander, S. K., Seinfeld, J. H., Enuiron. Sci. Technol., 3, 1175 Agency under grant number A P 1200-1, and in part at the Cali-
(1969). fornia Institute of Technology with the aid of a gift from the John
Hecht, T . A., Seinfeld, J. H., ibid., 6,47 (1972). A . McCarthy Foundation.

Volume 7, Number 3, March 1973 261