- Email: [email protected]
Theory of dispersion, transformation and deposition of atmospheric pollution using modified Green's functions
Vol. 13. pp. 239-246. Pergmon Press Ltd. 1979.Printed in Great Britain. ~rmospheric Enoironmenr
THEORY OF DISPERSION, TRANSFORMATION AND DEPOSITION OF ATMOSPHERIC POLLUTION USING MODIFIED GREEN’S FUNCTIONS G. ASTARITA
Chemical Engineering Department, University of Delaware Newark, DE 19711, U.S.A. J. WEI Chemical Engineering Department, Massachusetts Institute of Technology, Cambridge, MA 02193, U.S.A. and
G. IORIO* Istituto di Principi di Ingegneria Chimica, Universita di Napoli, Naples, Italy (First received 6 February 1978 and infinnlform
21 August 1978)
Abstract - Air pollution control strategies require a knowledge of the relation between the strength of the emission ofa primary pollutant at asource point and theconcentration ofa secondary or tertiary pollutant at a measure point. It is shown here that if the transformation and deposition kinetics follow first order dependence, the required relation can be given analytically as a tinite sum of modified Green’s functions. The general two-dimensional steady state case is given in detail. The approximations involved are physically well justified, and the computation of the spatial distribution of both primary and secondary pollutants can be performed on pocket calculators
INTRODUCTION
A primary air pollutant, such as sulfur dioxide, undergoes chemical transformation in the atmosphere to form secondary pollutants, such as sulfates and sulfuric acid. The pollutants may leave the atmosphere by wet or dry precipitation. The pollutants are also dispersed by the winds, and may be transported over hundreds of kilometers. Some acid rains in the northeast U.S. and in Scandinavia may have their origin in combustions in the industrial U.S. midwest and in England (Fisher, 1975; Likens, 1976). Consider a two-dimensional plane where there is a steady, spatially dispersed source of primary pollutant of strength S(x, y). After emission, the primary pollutant disperses by a diffusion mechanism with diffusivity D, is convected by wind with velocity u(x, y), and reacts to form a secondary pollutant with a first order kinetic rate constant k,. The equation of mass balance is (Reynolds et al., 1973) DV2Cl - u. VCr - k,Cr + S(x,y) = 0.
(1)
Similar equations, with more complex forms for the reaction term k,C1, can of course be written down for all the secondary pollutants. In general, the atmospheric pollution modeling problem can be formulated as follows: given an arbitrary (known) pollution strength distribution source S, and an arbitrary (known) wind speed distribution II, calculate the
concentration distributions C,(x, y) for both primary and secondary pollutants. Thanks to the linearity of the problem as described by Equation (l), a general closed form solution for any pollutant distribution source S(x,y) can be obtained, provided that a modified Green’s function for Equation (1) can beconstructed. In fact it is well known that the Green’s functions represent the solution of linear differential equations subjected to an impulse forcing function, i.e. for the case at hand, when S(x, y) = 6(x-x’, y-y’). Since the differential operator is linear, the response for an arbitrary foreign function can be then expressed by a superposition of impulses. Difficulties arise in connection with two issues. First, since arbitrary wind speed distributions u(x, y) need to be considered, the Green’s function to Equation (1) must be itself expressible as a functional of u(x, y). Second, if a Green’s function to the appropriate analogue of Equation (1) is to be constructed for secondary pollutants, arbitrary kinetic schemes need to be considered, and therefore the Green’s function will be a function of the specific kinetic scheme considered. Both these problems can be overcome, as is shown below, with reasonable assumptions regarding the types of wind speed distributions and chemical kinetics.
ONE-DIMENSIONAL STEADYSTATE SYSTEM
* Author to whom correspondence should be addressed.
We begin by considering only the primary pollutant 239
240
Ci ASTARITA.J.
WFI and G. IORIO
species in a one-dimensional steady state system. Equation (1) thus reduces to :
where both D and u are constants. If the source is a delta function S = &.X-X’), the solution is
ovAU-d 0
- %< .Y< .Y’ c = &!+> exp[IJ(x - 01 c = --‘---.exp[
400 distance,
km
Fig. t. Primary and secondary p~i~utants vs distance in onedimensional model. k, = Wh-‘. kz = 0.1 h‘~‘.
-y(x - X’)] X’
D(Il+ 1’) i
_
200
I
where
(4)
of course, crosswind diffusion cannot be neglected, since it is the only transport mechanism crosswind. This argument will be made precise later. We now turn out attention to the problem of secondary polutants. Consider a consecutive mechanism with first order kinetics
For a windless day, II = 0, the solutions reduce to p; (.Y- Y’) fr = 7 .~;~Texp \/ jj -t 1 I l/i, Cl = 2-&exp - J-g (x-u’) v
- c < .Y< \‘j (5) Y’<.x<
where C, is the concentration of the secondary pollutant, and k, describes the removal of the secondary pollutant by wet or dry deposition. Equation (2) is now given in matrix form as (Wei and Prater, 1962)
7
1
In the other limit of no dispersion, D = 0, the solutions reduce to [C,
-
=o
1 C1 = ;exp - k, -(Y--S’) 11 i
I
-
s’<.y<
.I
t6)
Since the matrix of rate constants can be diagonalized as
The R.H.S. of Equation (3) is the Green’s function to Equation (2), so that the general solution to Equation (2) with arbitrary S(.u’) is given by (Morse and Freshbach, 1953), c,(x)
=
.’i
G,(.+‘)S(x’)d.x’.
(7)
The form of the Green’s function G(.Y/s’) deserves some comment. For example, consider a set of conditions similar to SO1 in St. Louis; assume (Roth et u[., 1974; White et ul., 1976):
M
R
M -’
A._
and if we define (10) [tj
= ;,1,1
1 [I
u = 1.15km hh’ D=O.lgkm’h-’
then multiplying every term of Equation (8) from the left by M_’ yields
k, =O.lh-.’
one obtains p = 83.35 km-“, 7 = 0.0133 km-‘. The corresponding Green’s function G(.x/x’) is given in Fig. 1 (primary pollutant). The wind from the left overpowers the diffusion to. the left, so that the plume is essentially all on the right. This result suggests that, in the general case, diffusion in the downwind direction can be neglected;
D$[::]-u$j-[::
:;I fS
0
1 k,
=
(11)
. - __-._k,-k, 0 ! Ii] The solution to Cl is the same as before, but the solution to h2 is
241
Theory of dispersion, transformation and deposition of atmospheric pollution
The wind field is two-dimensional, and the independent variables to be considered are the Cartesian coordinates, x, y, of the region to be studied. Let u(x, y) be the wind velocity vector with components uX, uY.The equation of continuity is :
b2 =
b,
=
i
(12)
au
-2 ax
where
(13) y‘ =
&’
+-
au,
=o.
aY
(16)
If the viscosity of the air is neglected, the equation of motion take the form of the Euler equations, i.e. at steady state:
+ 4Dk, - u 20
’
Finally, C, is recovered from b, according to
where p is the pressure and p the density. In the absence of closed streamlines (i.e. when no large-scale whirls are present), the solution of Equations (16) and (17) is irrotational (Batchelor, 1970) i.e.
and
du, _ ay Pexp[fi’(x
- x’)]
au,
(18)
-K
In view of Equations (16) and (18) one may define a potential function $(x, y) and a stream function 4(x, y) such that: = G,(x, x’).
k, c2 = __ k, -k, +
-1
~ exp[ - y(x - x’)] D(P++Y)
i
l
1I
&expc-Y.(x - ~'11 , x’
J
J
(14)
Finally, for a spatially distributed source of primary pollutant S, the concentration of secondary pollutant is given by C,(x) =
G,(x, x’)S(x’)dx’.
(15)
For the same set of conditions discussed before, and with k, = 0.1 h-i (Rothetal., 1974; Whiteetal., 1976), one obtains $ = 83.34 km- ’ y’ = 0.0067 km- ‘. The corresponding Green’s function for the secondary pollutant, G,(xlx’), is given in Fig. 1. Again one may notice that the plume is essentially on the right. The procedure for transforming the coupled system of Equation (8) into the uncoupled system of Equation (11) can be generalized, as discussed below. Therefore, Green’s functions for secondary pollutants can be constructed for all chemical reaction systems that follow first order kinetics, reversible or irreversible.
TWO-DIMENSIONAL
IRROTATIONAL
WIND
Consider next a two-dimensional system, where vertical mixing is sufficiently rapid up to the inversion layer. A.E. 13;2- c
(19)
and Equations (16) and (18) are automatically satisfied. In Equation (19), U is a constant characteristic wind speed over the region considered, the introduction of which insures that the left-hand sides of Equation (19) are O(1). Therefore, also the right-hand sides are 0( 1). In the x, y plane, the curves 4 = constant are the wind streamlines, and the curves $ = constant are potential lines everywhere orthogonal to the streamlines. The normalization with U in Equation (19) ensures that A4 is of the order of the geometrical separation between two streamlines, and A$ is of the order of the geometrical separation between two potential lines. The x, y c* 4, $ transformation is everywhere smooth and invertible, since the Jacobian is non-zero unless the wind speed vanishes: det
wax way
u2+u; _ 1 =x__
I
[ a*laxa*iay
u2
q2'
(20)
The quantity q defined in Equation (20), is, ofcourse, a function of position ; its value is everywhere O(l), since it is the ratio of the squares of the actual and characteristic wind speeds. It will be seen that it is convenient to use the streamlines 4 = const. and the potential lines II/
G. ASTARITA,J. WEI and G. IORIO
242
= const. as coordinate lines. This allows easy mathematical formalization of the physical concepts of downwind (i.e. along a stream line) and crosswind (i.e. along a potential line). If the wind is rectilinear and accelerationless, the #, rj coordinate lines are Cartesian with the I// axis along the wind direction. In general, the wind may accelerate i.e. the stream lines need not be parallel, neither need they be rectilinear. It can be shown (Batchelor, 1970) that the crosswind derivative of the wind speed equals the ratio of the wind speed to the curvature radius of the streamline, R, (with the wind speed higher inside the bend). In terms of the #. $ coordinates and the function y, this result is expressed as follows
Equation (21) shows that, when the streamlines are rectilinear (with no need to be parallel) dq/@ = 0. It will be seen that the solution of the diffusion equations obtained below is rigorous in this case. In general, the derivative ~q/~~ is 0( l/&f. The rationale for introducing the considerations above is the following. If indeed under appropriate conditions the downwind diffusion can be neglected, the differential equations are best transformed from the x, y to the 4, $ coordinate system, since in the latter one coordinate (IL) varies only downwind, and the other one only crosswind. The appropriate differential equations are discussed in the next section. THE
DIFFERENTIAL
EQUATIONS
Let c = (et, . . . . c,) be the concentrations of the pollutants, which may be both primary or secondary ones. The rates of reaction r = (rl,. , rN) represent the rates of change of the concentrations in the absence of diffusion and convection ; they are in general functions of all the concentrations, say : r = r(c).
(22)
The differential equation for the steady state concentration distribution of pollutant i is:
pheric pollutants can be represented by Equation (24) is reasonable for chemical reactions involving a pollutant species (oxygen, light, rain) which is present at an excess or an essentially constant concentration. If Equation (24) is indeed applicable, we may follow the procedure introduced by Wei and Prater ( I Oh2 J.A new set of”equivalent concentrations” b = (h,.. h, I can be introduced, each element of which is a linear combination of the ci’s, and such that. in the absence of diffusion and convection I’S) where ki is an eigenvalue kinetic constant which can be calculated once the matrix K has been assigned. Since the differential Equations (24) are linear, and since the coefficients i?, I(, and uYappearing in them are the same for all, the equations can be transformed by the appropriate linear combination into the new set :
The differential Equations (26) are uncoupled and, therefore, can be solved independently. In the following, the index i will be dropped. We can now proceed to transform Equation (26) from the .Y,.I:coordinates to the 4, $ coordinates. The result is:
where it may be noticed that the vector-valued wind distribution u(x.3’) appears only through the scalarvalued wind speed distribution y(4,$); indeed, the vectorial character of the wind velocity is taken care of by the form of the I$, ri/ grid. It is now possible to proceed to an ordering argument which makes mathematically precise the concept of neglecting diffusion downwind. Let 15stand for the appropriate downwind length scale, and W for the appropriate crosswind length scale. If the concentration b is regarded as being O(l), which, since Equation (27) is linear, can always bJ; obtained by appropriate normalization, the four terms in Equation (27) have the following orders of magnitude. Diffusion crosswind
-ri((.,.
,c,).
(23)
The N differential equation such as Equation (23) are in genera1 coupled through the term r;, and their anaiytical solution is in the genera1 case a hopeless task. If the diffusivity is due to wind turbulence, it is necessarily the same for all pollutants. This fact allows the uncoupling of Equation (231, provided that the function t(c) is linear, say: r = x‘.c
Diffusion downwind (29) Convection (30)
(24)
with X a matrix of pseudo first-order kinetic constants. The assumption that the chemical kinetics of atmos-
Reaction rate (31)
Theory of dispersion, transformation and deposition of atmospheric pollution
Comparison of Equations mediately the result :
(30) and (31) yields im-
L = O(U/k)
(32)
i.e. the appropriate length scale downwind is the distance that the wind travels during the time of the chemical reaction l/k. Since the only mechanism of transport crosswind is diffusion, Equations (28) and (30) imply that :
w ---CO
(\i-! D
L
UL’
(33)
i.e. the ratio of the two length scales Wand L is of the order of the inverse square root of the Peclet number. Since with any reasonable combination of parameters values for U greater than 1 ms-i or so D/UL cc 1, the term DdZb/&,@, i.e. the downwind diffusion, can be neglected in Equation (27), thus reducing it to: d=b 8b Dw=U;i;L.+kqb.
OF EQUATION
Let g(&1/1) be the solution differential equation :
as can be verified by substitution into Equations (35) and (36). Now consider the following function :
h(#, $4 = s(rb,Wxp
The function h(e) satisfies the boundary condition :
h(h $‘I = s,(4)
of the following
Subject to the following boundary condition s(#, 9’) = s,(rb)
(36)
with g,(9) an arbitrary function. The solution g@, II/)is calculated as follows. First, a Green’s function G(& II/ I$‘, $‘) is obtained b y considering the case where g,(4) is a delta function :
a2h ah Dv-iI%-kqh=R,
= W-4’).
G = ,(~_~.)~PC-(d-~‘)‘~i4o(g_IL’)l
(37)
First of all one may notice that, if aq/a4 = 0 (i.e. the wind speed is constant crosswind), the residual R is zero. Therefore, the following conclusion is reached : when the wind streamlines are rectilinear (with no need to be parallel), the residual R is zero. Should the streamlines be curved, Equation (21) shows that the last two terms of the right-hand side of Equation (44) are O(D/Rz). In the following, we restrict attention to wind fields such that the radius of curvature of the streamhnes is not much smaller than the length scale L, say: R, 2 O(L).
R= -2~~exp[-~~~tqdJ]
* J
k a4 - -d$ VU JlL
X
- W?J’~
’ @Wt;. 5 --x
(39)
Since Equation (35) is linear, its solution under the general boundary condition (36) is :
(45)
Condition (45) is very realistic, since it simply implies that over the region considered the wind does not turn around 180”. If Equation (45) is satisfied the last two terms on the right-hand side of Equation (44) are at most O(D/Lz). Since the left-hand side of Equation (43) is O(U/L), the residual R can be taken, to within the same O(D/UL) for which Equation (34) holds, as:
(38)
as can easily be verified by substituting G(s) in Equations (35) and (37). The function If-) is the unit step function : l(t) =
(43)
with the residual R given by:
The Green’s function G(e) has the following form:
,/4M
(42)
and
(35)
Gt~,~~~‘,~‘)
(41)
(34)
Of course, since the diffusion downwind is of order D/UL with respect to the terms appearing in Equation (34), the solution of the latter needs to be sought only to within the same order of approximation.
SOLUTlON
243
+ 0(0/L’).
(46)
From Equations (38) and (40) one gets:
ag
d,
=
zz-
J
1 ac
_
I gd$W#
s zi
_
2
ac
abt’
= - Icgl/*=+~
h#Wd,
.#,‘=_*
+
Jx G$dg -X (47)
244
Ci. ASTARSA,
J. WEI and G.
IORIO
AN EXACT
and therefore : (4% Therefore, provided the kernel g,(@) in Equation (40) is sufficiently smooth that its derivative is at most 0(1/L), also the first term on the right-hand side of Equation (46) is 0(0/L’), and the entire residual R can be neglected to within the order of approximation D/VL. Hence, to within the same order /I(#+$) is a solution of Equation (34) which satisfies the boundary condition (42). To recapitulate, the residual R is exactly zero [i.e. h($, $) is an exact solution of Equation (34)) whenever the wind streamlines are rectilinear. If they are not rectilinear, under the additional restrictions imposed by Equation (45) and by the requirement that Idg,/dd, Inlax= 0(1/L), the residual R can be taken as zero to within the same order of approximation for which Equation (34) holds, i.e. provided the downwind diffusion is negligible. We are now in a position to solve the general problem where a distribution S(&+) of pollutant source densities is assigned (with S having units of mass of pollutant released per unit volume of air per unit time). Consider the differential surface between the potential lines $’ and $’ + dll/’ ; the amount of pollutant released over it is S(#,~‘)d~~, and this corresponds to a local increase of the pollutant’s concentration given by :
SOLliTlOh
In this section, we give actual calculations for a case where Equation (50) is an exact sotution of Equation (34). We consider the simplest chemical scheme which results in both a primary and a secondary pollutant, say : kt kl c, --+ c-2 --+CJ
(52)
with cr the primary and c2 the secondary pollutant : k I and k, are the actual kinetic constants, and entries of the matrix X in Equation (24). The transformation discussed following Equation (24) leads to the following equivalent concentration vector b = (b,,hz): h, = c,
(53)
h2 = (‘, + k?$c2 I
(54)
and the equivalent kinetic constants ki appearing in Equation (34) are: k, =k;
WI
k, = k;
06)
(notice that the identities (55) and (56) are due to the special simplicity of the chemical scheme considered). If H, and H2 are the modified Green’s function, as given by Equation ~5l)corr~ponding to b, and b,, the modified Green’s functions for the primary and secondary pollutant, H, and H,, are H,=
H,
(57)
H, = (H+I,);+~-. Correspondingly, the increase of pollutant concentration downwind (at $ > +‘) would be given by Equation (41), with db as given by Equation (49) substituted for g,(4). Since Equation (34) is linear, the total concentration of pollutants is simply the integral of all such differential concentrations. say :
where the function W(a) is:
158) 1
2
We consider a wind dist~bution with rectilinear (but not parallel) streamlines, so that the residual R in Equation (34) is exactly zero. In particular, with reference to Fig. 2, we consider a wind distribution which, downwind of the pollution source considered, can be approximated by a source flow. We refer to the point upwind where the downwind streamlines converge as the wind source. A cylindrical coordinate system Y,f?,centered at the wind source is chosen, and the pollution source is located at r = R, @ = 0. Of course, the functions r > R.
Hi need to be determined
only at
N&tir/#‘>$‘) Potbtion
The function H( * ) is a modified Green’s function for the probiem considered, which transforms the forcing function S(#, $‘) into the resulting effect b(#, $1; it is itself a function of the wind speed distribution function 4(& ti)*
Source
Fig. 2. Wind field and coordinate transformation from x, y to 4. CL.
Theory of dispersion, transformation and deposition of atmospheric pollution
245
From the known source flow velocity distribution, the 4, J/ grid is obtained as:
qQ=Rlni
(59)
4 = RB
(60)
where 0 = 0, $ = 0 have been assigned at the pol-
lution source. The function q is given by:
q=
ti
r
exp 0ii =0ii
where the characteristic velocity CJhas been taken as the wind speed at the pollution source. Equation (47) shows that q is constant along the potential lines and, therefore, R = 0. With straight forward algebra Equation (37) is reduced to:
Fig. 3. Primary pollutant concentration for wind fieldof. 2 and values in Table 1. Vertical scale: H = U,/4nD/K ; downwind scale: 1 grid interval = 1 km ; crosswind scale: 1 grid interval = 0.1 km.
(62) where : Bi = H,&i?iiii
(63)
P = 4D/UR
(64)
M, = &R/J.
(65)
Equations (57), (58) and (62) allow easy direct calculation of the primary and secondary pollutant distributions downwind of the pollution source. Actual calculations have been performed for the set of values of the relevant parameters reported in Table 1; the corresponding values of the dimensionless parameters P and Mi are given in the source table. Figures 3 and 4 are isometric plots of the normalized modified Green’s functions. When the distance R-+ x1, the source flow case degenerates to the case where the wind streamlines are both rectilinear and parallel; the 4, II/grid is in this case Cartesian, with the 4 = constant lines along the wind direction. Figures 5 and 6 are plots for this case, for the same group of values of the parameters as in Fig. 3, with of course R---t cc. It has to be noticed, for this particular case, that although values of physical parameters referring to a set of conditions similar to SO, in St. Louis have been adopted in the development of the example, no direct
Fig. 4. Secondary pollutant concentration ; same system as Fig. 3.
Table 1. Values of parameters U = lOkmh-’ D=0.25kmzh-’ R=lOkm k, = 1 h-’ k, =0.2h-’ P = 10-2 M, = 1 M, = 0.2 J4nURD = 17.7 km* h-’
Fig. 5. Primary pollutant concentration. R+ co, but all other values in Table 1.
246
G
ASTARITA.
J.
WFI
and G
IORIO
Work is under way to test the model performance effectiveness towards an actual example. presumably via a smoothing of the available experimental d:~ta
REFERE\\I(‘ES
Fig. 6. Secondary pollutant concentration : same system aS Fig. 5. comparison has been made. Indeed, data reported by Roth et a/. (1974) are so scattered that no matter how crude a model, comparison would be satisfactory provided the order of magnitude of the parameters is guessed correctly. Therefore a direct comparison between the experimental results obtained in the St. Louis study and the predictions ofthe model presented here cannot be performed In a significant way and merely qualitative conclusions can be drawn.
Batchelor G. K. (1970) An Inrroduc,tion to Fluid I>.vnumk.s Cambridge University Press. England. Fisher B. E. A. (1975) The long range transport oi sulfur dioxide. Atmospheric Enrironmem 9, 1063 1070. Likens G. E. (1976) Acid precipitation. Chem. Engnq Nurr’.\ November 22, 29 44. Morse P. M. and Freshbach A. (1953) Methods ofThrorrti~~~t/ Physics, Chaps. 7 and 12. McGraw-Hill, New York. Reynolds S. D., Roth P. M. and Seinfeld J. H. (1973) Mathematical modeling of photochemical air pollution I ; formulation of the model. .4tmospheric Encironmrnt 7, 1033 -1061. Roth P. M.. Roberts P. J. W.. Lin M. K.. Reynolds S. D. and Seinfeld J. H. (1974) Armosphrric~ Enoironmenf 8, 97 130. Wei J. and Prater C. D. (1962) The structure and analysis of complex reaction systems. Adr. (‘0tulysis 13, 203 397 White N. H.. Anderson J. A.. Blumenthal D. L.. Husar R. B.. Gillani N. V.. Husar J. D. and Wilson W. E. Jr.,(1976) Formation and transport of secondary air pollutants: ozone and aerosols in the St. Louis urban plume. Science 8, 187 198.
THEORY OF DISPERSION, TRANSFORMATION AND DEPOSITION OF ATMOSPHERIC POLLUTION USING MODIFIED GREEN’S FUNCTIONS G. ASTARITA
Chemical Engineering Department, University of Delaware Newark, DE 19711, U.S.A. J. WEI Chemical Engineering Department, Massachusetts Institute of Technology, Cambridge, MA 02193, U.S.A. and
G. IORIO* Istituto di Principi di Ingegneria Chimica, Universita di Napoli, Naples, Italy (First received 6 February 1978 and infinnlform
21 August 1978)
Abstract - Air pollution control strategies require a knowledge of the relation between the strength of the emission ofa primary pollutant at asource point and theconcentration ofa secondary or tertiary pollutant at a measure point. It is shown here that if the transformation and deposition kinetics follow first order dependence, the required relation can be given analytically as a tinite sum of modified Green’s functions. The general two-dimensional steady state case is given in detail. The approximations involved are physically well justified, and the computation of the spatial distribution of both primary and secondary pollutants can be performed on pocket calculators
INTRODUCTION
A primary air pollutant, such as sulfur dioxide, undergoes chemical transformation in the atmosphere to form secondary pollutants, such as sulfates and sulfuric acid. The pollutants may leave the atmosphere by wet or dry precipitation. The pollutants are also dispersed by the winds, and may be transported over hundreds of kilometers. Some acid rains in the northeast U.S. and in Scandinavia may have their origin in combustions in the industrial U.S. midwest and in England (Fisher, 1975; Likens, 1976). Consider a two-dimensional plane where there is a steady, spatially dispersed source of primary pollutant of strength S(x, y). After emission, the primary pollutant disperses by a diffusion mechanism with diffusivity D, is convected by wind with velocity u(x, y), and reacts to form a secondary pollutant with a first order kinetic rate constant k,. The equation of mass balance is (Reynolds et al., 1973) DV2Cl - u. VCr - k,Cr + S(x,y) = 0.
(1)
Similar equations, with more complex forms for the reaction term k,C1, can of course be written down for all the secondary pollutants. In general, the atmospheric pollution modeling problem can be formulated as follows: given an arbitrary (known) pollution strength distribution source S, and an arbitrary (known) wind speed distribution II, calculate the
concentration distributions C,(x, y) for both primary and secondary pollutants. Thanks to the linearity of the problem as described by Equation (l), a general closed form solution for any pollutant distribution source S(x,y) can be obtained, provided that a modified Green’s function for Equation (1) can beconstructed. In fact it is well known that the Green’s functions represent the solution of linear differential equations subjected to an impulse forcing function, i.e. for the case at hand, when S(x, y) = 6(x-x’, y-y’). Since the differential operator is linear, the response for an arbitrary foreign function can be then expressed by a superposition of impulses. Difficulties arise in connection with two issues. First, since arbitrary wind speed distributions u(x, y) need to be considered, the Green’s function to Equation (1) must be itself expressible as a functional of u(x, y). Second, if a Green’s function to the appropriate analogue of Equation (1) is to be constructed for secondary pollutants, arbitrary kinetic schemes need to be considered, and therefore the Green’s function will be a function of the specific kinetic scheme considered. Both these problems can be overcome, as is shown below, with reasonable assumptions regarding the types of wind speed distributions and chemical kinetics.
ONE-DIMENSIONAL STEADYSTATE SYSTEM
* Author to whom correspondence should be addressed.
We begin by considering only the primary pollutant 239
240
Ci ASTARITA.J.
WFI and G. IORIO
species in a one-dimensional steady state system. Equation (1) thus reduces to :
where both D and u are constants. If the source is a delta function S = &.X-X’), the solution is
ovAU-d 0
- %< .Y< .Y’ c = &!+> exp[IJ(x - 01 c = --‘---.exp[
400 distance,
km
Fig. t. Primary and secondary p~i~utants vs distance in onedimensional model. k, = Wh-‘. kz = 0.1 h‘~‘.
-y(x - X’)] X’
D(Il+ 1’) i
_
200
I
where
(4)
of course, crosswind diffusion cannot be neglected, since it is the only transport mechanism crosswind. This argument will be made precise later. We now turn out attention to the problem of secondary polutants. Consider a consecutive mechanism with first order kinetics
For a windless day, II = 0, the solutions reduce to p; (.Y- Y’) fr = 7 .~;~Texp \/ jj -t 1 I l/i, Cl = 2-&exp - J-g (x-u’) v
- c < .Y< \‘j (5) Y’<.x<
where C, is the concentration of the secondary pollutant, and k, describes the removal of the secondary pollutant by wet or dry deposition. Equation (2) is now given in matrix form as (Wei and Prater, 1962)
7
1
In the other limit of no dispersion, D = 0, the solutions reduce to [C,
-
=o
1 C1 = ;exp - k, -(Y--S’) 11 i
I
-
s’<.y<
.I
t6)
Since the matrix of rate constants can be diagonalized as
The R.H.S. of Equation (3) is the Green’s function to Equation (2), so that the general solution to Equation (2) with arbitrary S(.u’) is given by (Morse and Freshbach, 1953), c,(x)
=
.’i
G,(.+‘)S(x’)d.x’.
(7)
The form of the Green’s function G(.Y/s’) deserves some comment. For example, consider a set of conditions similar to SO1 in St. Louis; assume (Roth et u[., 1974; White et ul., 1976):
M
R
M -’
A._
and if we define (10) [tj
= ;,1,1
1 [I
u = 1.15km hh’ D=O.lgkm’h-’
then multiplying every term of Equation (8) from the left by M_’ yields
k, =O.lh-.’
one obtains p = 83.35 km-“, 7 = 0.0133 km-‘. The corresponding Green’s function G(.x/x’) is given in Fig. 1 (primary pollutant). The wind from the left overpowers the diffusion to. the left, so that the plume is essentially all on the right. This result suggests that, in the general case, diffusion in the downwind direction can be neglected;
D$[::]-u$j-[::
:;I fS
0
1 k,
=
(11)
. - __-._k,-k, 0 ! Ii] The solution to Cl is the same as before, but the solution to h2 is
241
Theory of dispersion, transformation and deposition of atmospheric pollution
The wind field is two-dimensional, and the independent variables to be considered are the Cartesian coordinates, x, y, of the region to be studied. Let u(x, y) be the wind velocity vector with components uX, uY.The equation of continuity is :
b2 =
b,
=
i
(12)
au
-2 ax
where
(13) y‘ =
&’
+-
au,
=o.
aY
(16)
If the viscosity of the air is neglected, the equation of motion take the form of the Euler equations, i.e. at steady state:
+ 4Dk, - u 20
’
Finally, C, is recovered from b, according to
where p is the pressure and p the density. In the absence of closed streamlines (i.e. when no large-scale whirls are present), the solution of Equations (16) and (17) is irrotational (Batchelor, 1970) i.e.
and
du, _ ay Pexp[fi’(x
- x’)]
au,
(18)
-K
In view of Equations (16) and (18) one may define a potential function $(x, y) and a stream function 4(x, y) such that: = G,(x, x’).
k, c2 = __ k, -k, +
-1
~ exp[ - y(x - x’)] D(P++Y)
i
l
1I
&expc-Y.(x - ~'11 , x’
J
J
(14)
Finally, for a spatially distributed source of primary pollutant S, the concentration of secondary pollutant is given by C,(x) =
G,(x, x’)S(x’)dx’.
(15)
For the same set of conditions discussed before, and with k, = 0.1 h-i (Rothetal., 1974; Whiteetal., 1976), one obtains $ = 83.34 km- ’ y’ = 0.0067 km- ‘. The corresponding Green’s function for the secondary pollutant, G,(xlx’), is given in Fig. 1. Again one may notice that the plume is essentially on the right. The procedure for transforming the coupled system of Equation (8) into the uncoupled system of Equation (11) can be generalized, as discussed below. Therefore, Green’s functions for secondary pollutants can be constructed for all chemical reaction systems that follow first order kinetics, reversible or irreversible.
TWO-DIMENSIONAL
IRROTATIONAL
WIND
Consider next a two-dimensional system, where vertical mixing is sufficiently rapid up to the inversion layer. A.E. 13;2- c
(19)
and Equations (16) and (18) are automatically satisfied. In Equation (19), U is a constant characteristic wind speed over the region considered, the introduction of which insures that the left-hand sides of Equation (19) are O(1). Therefore, also the right-hand sides are 0( 1). In the x, y plane, the curves 4 = constant are the wind streamlines, and the curves $ = constant are potential lines everywhere orthogonal to the streamlines. The normalization with U in Equation (19) ensures that A4 is of the order of the geometrical separation between two streamlines, and A$ is of the order of the geometrical separation between two potential lines. The x, y c* 4, $ transformation is everywhere smooth and invertible, since the Jacobian is non-zero unless the wind speed vanishes: det
wax way
u2+u; _ 1 =x__
I
[ a*laxa*iay
u2
q2'
(20)
The quantity q defined in Equation (20), is, ofcourse, a function of position ; its value is everywhere O(l), since it is the ratio of the squares of the actual and characteristic wind speeds. It will be seen that it is convenient to use the streamlines 4 = const. and the potential lines II/
G. ASTARITA,J. WEI and G. IORIO
242
= const. as coordinate lines. This allows easy mathematical formalization of the physical concepts of downwind (i.e. along a stream line) and crosswind (i.e. along a potential line). If the wind is rectilinear and accelerationless, the #, rj coordinate lines are Cartesian with the I// axis along the wind direction. In general, the wind may accelerate i.e. the stream lines need not be parallel, neither need they be rectilinear. It can be shown (Batchelor, 1970) that the crosswind derivative of the wind speed equals the ratio of the wind speed to the curvature radius of the streamline, R, (with the wind speed higher inside the bend). In terms of the #. $ coordinates and the function y, this result is expressed as follows
Equation (21) shows that, when the streamlines are rectilinear (with no need to be parallel) dq/@ = 0. It will be seen that the solution of the diffusion equations obtained below is rigorous in this case. In general, the derivative ~q/~~ is 0( l/&f. The rationale for introducing the considerations above is the following. If indeed under appropriate conditions the downwind diffusion can be neglected, the differential equations are best transformed from the x, y to the 4, $ coordinate system, since in the latter one coordinate (IL) varies only downwind, and the other one only crosswind. The appropriate differential equations are discussed in the next section. THE
DIFFERENTIAL
EQUATIONS
Let c = (et, . . . . c,) be the concentrations of the pollutants, which may be both primary or secondary ones. The rates of reaction r = (rl,. , rN) represent the rates of change of the concentrations in the absence of diffusion and convection ; they are in general functions of all the concentrations, say : r = r(c).
(22)
The differential equation for the steady state concentration distribution of pollutant i is:
pheric pollutants can be represented by Equation (24) is reasonable for chemical reactions involving a pollutant species (oxygen, light, rain) which is present at an excess or an essentially constant concentration. If Equation (24) is indeed applicable, we may follow the procedure introduced by Wei and Prater ( I Oh2 J.A new set of”equivalent concentrations” b = (h,.. h, I can be introduced, each element of which is a linear combination of the ci’s, and such that. in the absence of diffusion and convection I’S) where ki is an eigenvalue kinetic constant which can be calculated once the matrix K has been assigned. Since the differential Equations (24) are linear, and since the coefficients i?, I(, and uYappearing in them are the same for all, the equations can be transformed by the appropriate linear combination into the new set :
The differential Equations (26) are uncoupled and, therefore, can be solved independently. In the following, the index i will be dropped. We can now proceed to transform Equation (26) from the .Y,.I:coordinates to the 4, $ coordinates. The result is:
where it may be noticed that the vector-valued wind distribution u(x.3’) appears only through the scalarvalued wind speed distribution y(4,$); indeed, the vectorial character of the wind velocity is taken care of by the form of the I$, ri/ grid. It is now possible to proceed to an ordering argument which makes mathematically precise the concept of neglecting diffusion downwind. Let 15stand for the appropriate downwind length scale, and W for the appropriate crosswind length scale. If the concentration b is regarded as being O(l), which, since Equation (27) is linear, can always bJ; obtained by appropriate normalization, the four terms in Equation (27) have the following orders of magnitude. Diffusion crosswind
-ri((.,.
,c,).
(23)
The N differential equation such as Equation (23) are in genera1 coupled through the term r;, and their anaiytical solution is in the genera1 case a hopeless task. If the diffusivity is due to wind turbulence, it is necessarily the same for all pollutants. This fact allows the uncoupling of Equation (231, provided that the function t(c) is linear, say: r = x‘.c
Diffusion downwind (29) Convection (30)
(24)
with X a matrix of pseudo first-order kinetic constants. The assumption that the chemical kinetics of atmos-
Reaction rate (31)
Theory of dispersion, transformation and deposition of atmospheric pollution
Comparison of Equations mediately the result :
(30) and (31) yields im-
L = O(U/k)
(32)
i.e. the appropriate length scale downwind is the distance that the wind travels during the time of the chemical reaction l/k. Since the only mechanism of transport crosswind is diffusion, Equations (28) and (30) imply that :
w ---CO
(\i-! D
L
UL’
(33)
i.e. the ratio of the two length scales Wand L is of the order of the inverse square root of the Peclet number. Since with any reasonable combination of parameters values for U greater than 1 ms-i or so D/UL cc 1, the term DdZb/&,@, i.e. the downwind diffusion, can be neglected in Equation (27), thus reducing it to: d=b 8b Dw=U;i;L.+kqb.
OF EQUATION
Let g(&1/1) be the solution differential equation :
as can be verified by substitution into Equations (35) and (36). Now consider the following function :
h(#, $4 = s(rb,Wxp
The function h(e) satisfies the boundary condition :
h(h $‘I = s,(4)
of the following
Subject to the following boundary condition s(#, 9’) = s,(rb)
(36)
with g,(9) an arbitrary function. The solution g@, II/)is calculated as follows. First, a Green’s function G(& II/ I$‘, $‘) is obtained b y considering the case where g,(4) is a delta function :
a2h ah Dv-iI%-kqh=R,
= W-4’).
G = ,(~_~.)~PC-(d-~‘)‘~i4o(g_IL’)l
(37)
First of all one may notice that, if aq/a4 = 0 (i.e. the wind speed is constant crosswind), the residual R is zero. Therefore, the following conclusion is reached : when the wind streamlines are rectilinear (with no need to be parallel), the residual R is zero. Should the streamlines be curved, Equation (21) shows that the last two terms of the right-hand side of Equation (44) are O(D/Rz). In the following, we restrict attention to wind fields such that the radius of curvature of the streamhnes is not much smaller than the length scale L, say: R, 2 O(L).
R= -2~~exp[-~~~tqdJ]
* J
k a4 - -d$ VU JlL
X
- W?J’~
’ @Wt;. 5 --x
(39)
Since Equation (35) is linear, its solution under the general boundary condition (36) is :
(45)
Condition (45) is very realistic, since it simply implies that over the region considered the wind does not turn around 180”. If Equation (45) is satisfied the last two terms on the right-hand side of Equation (44) are at most O(D/Lz). Since the left-hand side of Equation (43) is O(U/L), the residual R can be taken, to within the same O(D/UL) for which Equation (34) holds, as:
(38)
as can easily be verified by substituting G(s) in Equations (35) and (37). The function If-) is the unit step function : l(t) =
(43)
with the residual R given by:
The Green’s function G(e) has the following form:
,/4M
(42)
and
(35)
Gt~,~~~‘,~‘)
(41)
(34)
Of course, since the diffusion downwind is of order D/UL with respect to the terms appearing in Equation (34), the solution of the latter needs to be sought only to within the same order of approximation.
SOLUTlON
243
+ 0(0/L’).
(46)
From Equations (38) and (40) one gets:
ag
d,
=
zz-
J
1 ac
_
I gd$W#
s zi
_
2
ac
abt’
= - Icgl/*=+~
h#Wd,
.#,‘=_*
+
Jx G$dg -X (47)
244
Ci. ASTARSA,
J. WEI and G.
IORIO
AN EXACT
and therefore : (4% Therefore, provided the kernel g,(@) in Equation (40) is sufficiently smooth that its derivative is at most 0(1/L), also the first term on the right-hand side of Equation (46) is 0(0/L’), and the entire residual R can be neglected to within the order of approximation D/VL. Hence, to within the same order /I(#+$) is a solution of Equation (34) which satisfies the boundary condition (42). To recapitulate, the residual R is exactly zero [i.e. h($, $) is an exact solution of Equation (34)) whenever the wind streamlines are rectilinear. If they are not rectilinear, under the additional restrictions imposed by Equation (45) and by the requirement that Idg,/dd, Inlax= 0(1/L), the residual R can be taken as zero to within the same order of approximation for which Equation (34) holds, i.e. provided the downwind diffusion is negligible. We are now in a position to solve the general problem where a distribution S(&+) of pollutant source densities is assigned (with S having units of mass of pollutant released per unit volume of air per unit time). Consider the differential surface between the potential lines $’ and $’ + dll/’ ; the amount of pollutant released over it is S(#,~‘)d~~, and this corresponds to a local increase of the pollutant’s concentration given by :
SOLliTlOh
In this section, we give actual calculations for a case where Equation (50) is an exact sotution of Equation (34). We consider the simplest chemical scheme which results in both a primary and a secondary pollutant, say : kt kl c, --+ c-2 --+CJ
(52)
with cr the primary and c2 the secondary pollutant : k I and k, are the actual kinetic constants, and entries of the matrix X in Equation (24). The transformation discussed following Equation (24) leads to the following equivalent concentration vector b = (b,,hz): h, = c,
(53)
h2 = (‘, + k?$c2 I
(54)
and the equivalent kinetic constants ki appearing in Equation (34) are: k, =k;
WI
k, = k;
06)
(notice that the identities (55) and (56) are due to the special simplicity of the chemical scheme considered). If H, and H2 are the modified Green’s function, as given by Equation ~5l)corr~ponding to b, and b,, the modified Green’s functions for the primary and secondary pollutant, H, and H,, are H,=
H,
(57)
H, = (H+I,);+~-. Correspondingly, the increase of pollutant concentration downwind (at $ > +‘) would be given by Equation (41), with db as given by Equation (49) substituted for g,(4). Since Equation (34) is linear, the total concentration of pollutants is simply the integral of all such differential concentrations. say :
where the function W(a) is:
158) 1
2
We consider a wind dist~bution with rectilinear (but not parallel) streamlines, so that the residual R in Equation (34) is exactly zero. In particular, with reference to Fig. 2, we consider a wind distribution which, downwind of the pollution source considered, can be approximated by a source flow. We refer to the point upwind where the downwind streamlines converge as the wind source. A cylindrical coordinate system Y,f?,centered at the wind source is chosen, and the pollution source is located at r = R, @ = 0. Of course, the functions r > R.
Hi need to be determined
only at
N&tir/#‘>$‘) Potbtion
The function H( * ) is a modified Green’s function for the probiem considered, which transforms the forcing function S(#, $‘) into the resulting effect b(#, $1; it is itself a function of the wind speed distribution function 4(& ti)*
Source
Fig. 2. Wind field and coordinate transformation from x, y to 4. CL.
Theory of dispersion, transformation and deposition of atmospheric pollution
245
From the known source flow velocity distribution, the 4, J/ grid is obtained as:
qQ=Rlni
(59)
4 = RB
(60)
where 0 = 0, $ = 0 have been assigned at the pol-
lution source. The function q is given by:
q=
ti
r
exp 0ii =0ii
where the characteristic velocity CJhas been taken as the wind speed at the pollution source. Equation (47) shows that q is constant along the potential lines and, therefore, R = 0. With straight forward algebra Equation (37) is reduced to:
Fig. 3. Primary pollutant concentration for wind fieldof. 2 and values in Table 1. Vertical scale: H = U,/4nD/K ; downwind scale: 1 grid interval = 1 km ; crosswind scale: 1 grid interval = 0.1 km.
(62) where : Bi = H,&i?iiii
(63)
P = 4D/UR
(64)
M, = &R/J.
(65)
Equations (57), (58) and (62) allow easy direct calculation of the primary and secondary pollutant distributions downwind of the pollution source. Actual calculations have been performed for the set of values of the relevant parameters reported in Table 1; the corresponding values of the dimensionless parameters P and Mi are given in the source table. Figures 3 and 4 are isometric plots of the normalized modified Green’s functions. When the distance R-+ x1, the source flow case degenerates to the case where the wind streamlines are both rectilinear and parallel; the 4, II/grid is in this case Cartesian, with the 4 = constant lines along the wind direction. Figures 5 and 6 are plots for this case, for the same group of values of the parameters as in Fig. 3, with of course R---t cc. It has to be noticed, for this particular case, that although values of physical parameters referring to a set of conditions similar to SO, in St. Louis have been adopted in the development of the example, no direct
Fig. 4. Secondary pollutant concentration ; same system as Fig. 3.
Table 1. Values of parameters U = lOkmh-’ D=0.25kmzh-’ R=lOkm k, = 1 h-’ k, =0.2h-’ P = 10-2 M, = 1 M, = 0.2 J4nURD = 17.7 km* h-’
Fig. 5. Primary pollutant concentration. R+ co, but all other values in Table 1.
246
G
ASTARITA.
J.
WFI
and G
IORIO
Work is under way to test the model performance effectiveness towards an actual example. presumably via a smoothing of the available experimental d:~ta
REFERE\\I(‘ES
Fig. 6. Secondary pollutant concentration : same system aS Fig. 5. comparison has been made. Indeed, data reported by Roth et a/. (1974) are so scattered that no matter how crude a model, comparison would be satisfactory provided the order of magnitude of the parameters is guessed correctly. Therefore a direct comparison between the experimental results obtained in the St. Louis study and the predictions ofthe model presented here cannot be performed In a significant way and merely qualitative conclusions can be drawn.
Batchelor G. K. (1970) An Inrroduc,tion to Fluid I>.vnumk.s Cambridge University Press. England. Fisher B. E. A. (1975) The long range transport oi sulfur dioxide. Atmospheric Enrironmem 9, 1063 1070. Likens G. E. (1976) Acid precipitation. Chem. Engnq Nurr’.\ November 22, 29 44. Morse P. M. and Freshbach A. (1953) Methods ofThrorrti~~~t/ Physics, Chaps. 7 and 12. McGraw-Hill, New York. Reynolds S. D., Roth P. M. and Seinfeld J. H. (1973) Mathematical modeling of photochemical air pollution I ; formulation of the model. .4tmospheric Encironmrnt 7, 1033 -1061. Roth P. M.. Roberts P. J. W.. Lin M. K.. Reynolds S. D. and Seinfeld J. H. (1974) Armosphrric~ Enoironmenf 8, 97 130. Wei J. and Prater C. D. (1962) The structure and analysis of complex reaction systems. Adr. (‘0tulysis 13, 203 397 White N. H.. Anderson J. A.. Blumenthal D. L.. Husar R. B.. Gillani N. V.. Husar J. D. and Wilson W. E. Jr.,(1976) Formation and transport of secondary air pollutants: ozone and aerosols in the St. Louis urban plume. Science 8, 187 198.