Numerical simulation of air pollution in urban areas: Model development

ooo4-6981/a7 53.00+0.00 Fwpmm Joumlc Ltd.

~nmwp4wricEwirmmmr Vol. 21. No. I. pp. 201-212. 1987. Fnnted in Goal Britain.

NUMERICAL

SIMULATION OF AIR POLLUTION AREAS: MODEL DEVELOPMENT

IN URBAN

JIA-YEONGKu and S. TRIWKRAMA RAO Division of Air Resources,New York State Department of Environmental Conservation. Albany, NY 12233, U.S.A.

and K. SHANRAR RAO Atmospheric Turbulena and Diffusion Division, National Daanic and Atmospheric Administration, Oak Ridge, TN 37830, U.S.A. (First receiwd 20 January 1986 cm&infinof/orm 7 July 1986) Abstract-A threedimensional, grid-based numerical air pollution model for the estimation of air pollutant concentrations in an urban area is developed. Based on the continuity equation, the modeling system incorporates the combined influenas of advcctivc transport, turbulent diffusion, chemical transformation, soura emissionsand surface removal of air contaminants. Recent developments in plume rise and plume penetration processes, objective wind field analysis procedures and numerical solution techniques incorporated into the model are described. Key word index: Urban air pollution models, numerical grid model, objective wind field analysis, numerical techniques, turbulent diffusion. deposition, chemical transformation.

I. INTRODUCTION Atmospheric diffusion models are being widely used to study the complicated relationships between air quality and emission sources as a function of meteorological and source conditions. Mathematical models used in the air pollution field range from simple empirical models to very complex numerical models. In general, these models are based on the equation governing the pollutant concentration consistent with the physical principle of mass conservation. The governing equation is solved with the initial and boundary conditions for a given source emission rate. However, due to the complexity of the atmospheric circulation patterns in the urban areas and inadequate data availabihty, various empirical approaches have been sought for practical purposes, based on either statistical theory, or observed features, or sometimes even arbitrary assumptions. Most of the currently available models utilize the Turner (1964) criteria for stability classification. These produce a bias toward neutral atmospheric conditions, which affects the estimation of plume rise and plume dispersion. Recently, substantial progress has been made tow;trd the understanding of the proccsscs governing the mean and turbulence structure of the planetary boundary layer (PBL). The PBL theory was improved for the convective conditions through numerical modeling (Deardorff, 1972). field observations (Kaimal et 01.. 1976), and laboratory experiments (Willis and htrdorll; 1978). Shir and Shich (1974) 201

developed a three-dimensional grid model for the urban environment and utilized an objective analysis scheme for the preparation of input data. The model was tested with SOz data from the St Louis metropolitan area and the model predictions agreed favorably with experimental measurements for both long term and short-term concentrations. Weil and Brower (1984) demonstrated that a Gaussian model with upto-date PBL theory can provide better results. The object of this paper is to develop a numerical air pollution model incorporating the recent advances in PBL dynamics, with emphasis on the daytime convective conditions which usually yield the highest ground level concentration from elevated buoyant sources. The model incorporates: (1) bulk Richardson number as the stability parameter instead of the Turner criteria to determine the atmospheric stability; (2) an eddy diffusivity profile for the convective boundary layer based on a numerical turbulence model (Lamb and Durran, 1978); (3) Briggs’ (1975, 1984) new plume rise formulae for convective conditions, and new methods for plume penetration into elevated stable layers as opposed to the usual “all or none” type of treatment. In addition, since wind field interpolation schemes and numcricttl tnctltods phty itn important role in dctermining the pollutant transport, a nondivergent wind field interpolation schemewhich produces a wind field with the structure inherent in the original wind measurements is included in the model in order to better represent the urban wind field charaeteristies and

IO ttutintnin

the

physical

constraint

of

mass

202

JIA-YEONG

conservation. The model also utilizes the method of fractional steps for the numerical solution of the turbulent diffusion equation. The advection terms are treated with the flux-limited method to minimize the artificial loss or gain of the total mass within the modeling domain. Further, the surface removal of the pollutants is considered by incorporating a deposition velocity in the surface boundary condition. A linear chemical transformation rate process is also included in the numerical model for the treatment of chemical sources/sinks.

Ku er al

layer. The height H is always greater than or equal to the mixed layer depth in the model formulation to allow plumes trapped in the stable layer to beentrained into the mixed layer during the late morning hours when the mixed layer rises rapidly. The initial conditions are arbitrarily set to zero, because it is difficult to determine the initial conditions at all grid points due to the sparsity of data. The initial conditions have been found to be important only in the initial period of simulation (Shir and Shieh, 1974; McRae et al., 1982).

2. MODEL FORMULATION

2.3. Treatment oj meteorological variables 2.1. Governing equation The mathematical basis for the urban air pollution

model is the conservation of mass equation, which describes advection, turbuknt diflitsion,chemical reaction, removal and emissions of the pollutants. In this study, a first order chemical transformation from sulfur dioxide (SOJ to particulate sulfate (SO: - ) is considered. The ensemble-averaged quation for the concentration can be expressed as follows:

i = 1 and 2

(2.1)

where C, and C2 are concentrations of SOz and sulfate, respectively; t is time; u, 1)and ware mean wind components along the x, y and z coordinates; K, and K, are horizontal and vertical eddy diffusion coefficients; Ri and Qi are terms characterizing the chemical source/sink and direct emission sources for species i, respectively. The chemical transformation from SO2 to sulfate can be expressed by R, = -k,C, and Rz = yk.C,, where k,, is reaction rate constant and y is the ratio of molecular weight of sulfate to SOI.

The variables required for the integration of Equation (2.1) and its boundary conditions are wind field, eddy diffusivities, chemical transformation rate, deposition velocities, emission rates and mixing height. The methods used to obtain the values of input variables for the model are discussed below. 2.3.1. The windfield. A mass-consistent wind field is essential in the air quality simulation model to eliminate artificial loss or generation in the species mass balance. A divergence-free objective analysis scheme for the St Louis urban area developed by Clark and Eskridge (1977) has been applied in preparing the wind fields for the model simulation. This scheme produces a wind field with the structure inherent in the original wind measurements and maintains the physical constraint of mass conservation. Two steps are involved in the numerical scheme. The first step is to interpolate the measuring stations data to the grid points, The weighting is a function of the distance between the station and the grid point. Xi,j=Xi,j+‘;;;;m,k=l,2

N (2.3)

where N is the number of data points within the scan circlecentered at the grid point (i,j); W f+jisa weighting factor and d’*“’ is a correction factor determined at the mth iteration of the station point. The weighting function, W f, j (r) is defined as,

[(Rf- (ri.jY)/(Rf + (ri.j)2)]4v r;,j

G R,

(2

4)

0, rt., > R,

2.2. Boundary conditions The boundary conditions are, - K.dCJ&

,...,

(8J

= Fi(x, y, t) at z = 0

Z0

atz=H

K&fax

= 0

at x = xW or xr

K,aCi/dy

= 0

aty=ysoryN

(2.2)

where F, = v& (0) are the surface deposition fluxes, and u,{ is the dry deposition velocity of pollutant i, xw and xE are the west and east, and ys and yN are south and north boundaries for the modeling domain. The height of tbe top boundary, H, is varied in accordance with the growth of the depth of the convective mixed

where R, is the scan radius, and r is the distance from grid point (i, j) to station k. The value of the correction factor (d”, “) is determined by subtracting the observed data value from the previous iteration value at a point corresponding to the location of the station. The second step is the adjustments of u and v components to reduce the divergence or convergence in the wind field. The adjustment velocity is estimated based on the divergence calculated at (i, j): Di.j=

(Ui+l,j-Ui-I,j)/(Xi*l

-xI-L) (2.5)

+(ui,j+I-ui,j-l)/(Yj+l-Yj-1).

Liu

and

Goodin

(1976) defined

the

following

re~tion~i~~

changeswith ~tmosp~r~ stabihty and surface rougb-

s+ t

Qf+t,j+A*l,J4,j

ui+l,j" g+’ I

1.1

t+ I “i.j+l 1+1 Ui,j_1

= uf-I.J*A-l,J$.J = u;,j* 1 +$,j* =Uf,jd

1

-12.

I ?5,j $

c

1.

6,

(2.61

J

where 5%j is a weighting factor, ul end u*are the sth iteration V&ES, and r7i.i and &,j are adjustment velocities The empiricat values for& from Chsstkand E&ridge (1977)are: (a) the grid point nearest to the ~~~to~ statiun is treitted as center point and assigned a weight of 035; (b) the ad&ent grid po%ts are assigned a weight of 0.5 and utbcr points are assigned a weight of 1.0, Substituting (2.6) into (2s) and partitioning Di,j into the u and u components based on the ratio of the weighting factor, the adjustment velocity components can be derived as, c Ui.j

=

-Di.j(xi+l--xi,,)l(fr+*.j+~-I,j +fi.j+

Ui_j=

1 +_fi.j-1)

-Di,j(Yj+r-~jYj,if/tf;lcI,j-f-fi-t,j +fi,j*X

+fi,

ness acoordittg to the v&es given in Table 1 (Irwin, 1979). 2.32, Atmsphsric stability. The rates of diffusion of the plume in the atmosphere are highly dependent on the atmospheric stability. Usually the stability classification proposed by Pasquill (1962),which depends on the solar insolation, wind speed and cloud cover, is employed in many of the currently used air pollution models. WeiIand Prower (1984)argued that the method was too subjective and yields too many neutral ~~i~~~ fn t&s secrion, tftedisassiort is focused on the use of the surface stability ~rnet~~ and the detc~~tjon of Moni~bukhov Itmgth (t) and friction velocity (u,), to estimate the eddy diffusivities. Direct determination of the t from its definition requires measurements of heat flux and friction velocity. Unfortunately, such measurements are not routinely available.Therefore, Land u* are derived by

(2.71

j- I)*

~~t~~~~ @S]s {2.&fand (2.7) awe tkn sobed it+z~~ti&y, untif the maximum divergence on the grid is

equal to zero or reduced to an acceptable level (lo-‘s-‘, as in &here, 1981). Figure 1 shows an example of maximum divergent as a function of iteration count. Upper level winds are obtained by extrapolation of the surface wind field using the power law profile, Ivl = i%if%/z,)”

(2.81

where rr, is the wind vector at s,, and the exponent p

Roughness (1.01m

stability

P

O.iOm P

1.00m P

3L?Q m P

A

o.as

0.21

It::

0.08 O.ll 0.09

0.17

:

0.20 0.17

0.31 0.28

; F

0.12 0,34 O.f3

0.16 0.32 0.54

0.27 0.38 0.6I

0.37 0.47 0.69

JIA-YEONG Ku er al.

204

using the bulk Richardson number, and expressions for the velocity and temperature profiles in the surface layer. The bulk Richardson number is defined by (2.9) Theatmosphere is stable with positive &and unstable with negative Rib. The velocity and temperature profiles within the surface layer are based on the expressions given by Businger et al. (1971). 1 + 4.7 i/L, neutral and stable (l-15

z/L)-‘/4,unstable

(2.10)

0.74 + 4.7 z/L, neutral and stable 0.74 (I - 9 r/L)- ‘!z, unstable. (2.1 I) For neutral and stable conditions, an analytic solution for L by solving Equations (2.9). (2.10) and (2.11) is,

AZ/L = ln(zJz,)

velocity scale, w,. The latter is given by (2.16)

A value ofk = 0.35 for the von Kirman constant given by Businger et al. (1971) is used in this study. 2.3.3. Eddydifrusiviries. Central to the treatment of the turbulent ditTusion process in Equation (2.1) is the estimation of the horizontal and vertical diffusivity coefficients, K, and K,. The parameterization of the turbulent flux from K-theory usually is based on mixing length theory. This theory was modified by dimensional arguments and empirical stability considerations to represent the planetary boundary layer physics (Blackadar, 1962; Shir, 1973; Yu, 1977). K,. is a function of vertical wind shear and lapse rate, but the functional relationship between K, and these variables is largely unknown. A method which employs only

9.4 Rib - 0.74 + (4.8888 R,b + 0.5476) ’ ‘* 9.4-44.18Rib

For the unstable condition, AZ/L is a transcendental function. A least-square-fitting method to compute L (Chang, 1978) is, Y = a,+a,x+a2x’+a3x3

(2.13)

where Y = In ( - z2/L); x = In ( - R,,,,) and Ri, is the modified bulk Richardson number, defined as Rim = Rib (zz/Az).

(2.14)

Here, z1 and z2 are the heights above the surface within the surface boundary layer; AZ, Au and A0 are the differences in height, wind speed and temperature between z2 and z,, respectively. The coefficients ao, a,, al and a3 are functions of :Jz, only:

10-3q-0.10311

a2 = 0.59585 x lo-‘$0.12886

x

(2.12)

those available parameters (e.g., L, surface roughness Zc, and u*) to determine the diffusion coefficients is utilized in this study. (a) Vertical eddy diffusivity for unstable conditions. Vertical transport of pollutant material is dominated by turbulent diffusion under convective conditions. Yamada (1977) has found diffusivities on the order of 100 m* s-’ in his numerical modeling studies. It is generally known that K, in the CBL increases with height at a rate somewhat greater than linear in the lower levels. The region of largest K, from 0.2 < z/Z, < 0.55 occursat a height where X’/Sz is nearly zero. At

10-2q3

u0 =0.6351+0.93823q-0.13151q2+0.73386x a, =0.97538-0.21956x

.

x lo-*q*+O.81699~

lo-‘q-O.36299

10e4q3

x lo-‘q* -0.3422 x 10W5q3

cj = 0.9037 x 10 - 4 + 0.26283 x 10 -‘q + 0.49967 x 10 - ‘q* - 0.76682 x 10 - ‘q3 and 9 = In (zJz,).

The value for u, is calcukrtcd by intcgmting Equation (2. IO) 3% u, = kuk)/J/,,

(2.15)

where expressions for the integral, I/I,,,are summarized in Table 2. For the convcctivc boundary layer, the two most important variables conirolling turbulence are convective boundary layer (CBL) depth, Z, and convective

the top of the mixed layer, K, approaches zero (Crane et al., 1977; Kaimal er al.. 1976). The vcrti~ldill\rsivily ofpollul;u+i is though1 lo bc more strongly related to the eddy diflusivity of heat rather than momentum (Crane et a/., 1977). In the

surface layer (0 6 z/Z; < 0.05). it is rccognizcd thal h’,. is a function of z/L only. The expression suggested by McRae er ul. (1982) is used in the numerical model for

205

Numerical simulation of air pollution in urban areas: model development

Table 2. Momentum integrals +, (z/L) for different stability conditions $, = jkdzlz Stable

ln(z/zO)+4.7(z-2,)/L

Neutral

In (z/to)

Unstable In

0 G Z/Zi

C

K”

(l-ISz/L)~+I (I-lSz/Ljf-1

0.05: = 2.5 w,Zj(kz/Z;)4’3

(1 - 15 z/L)*.

(2.17)

Above the surface layer, there is considerable uncertainty regarding the magnitude and height dependence of the eddy diffusivity. Lamb and Durran (1978) derived an empirical expression for K, using the threedimensional numerical turbulence mode1 of Deardorff (1972). When scaled with w, and Zi, the diffusivity profile can be expressedas a function of (z/Z& K,/‘w*Z, = 0.021+ O.~8(z/Zi)+

1*352(z/Zi)’

-4.096(z/‘Zi)3 + 2.56(z/Zif4.

(2.18)

According to (2.i8), the maximum value of K, occurs at z/Z; 1: 0.5 and has a magnitude (0.21 w,ZJ. For typical convective conditions, w* = 1.5 m s - ’ and Zi = 1000 m, this corresponds to a diffusivity of order of

magnitude 300 m2 s- ‘. Above the maximum point, K,, decreasesto a small value near Z, where turbulence transport is very small. There are some limitations in applying Ic, in convection conditions. Due to small values of K, at upper levels, the downwind displacement of pollutants from elevated sources is not handled properly. This may be compensated by reducing the plume rise or decreasing the physical stack height. For a ground level source, K,,, diffusivity works well up to a distance of x % 0.5 (u/w*)Zi, It will overpredict the rapidly declining surface concentrations beyond this point [from X * 0.5 (u/w,)Z~ to X x 2.5(U/w,)Zi], due to a

temporary upward flux of pollutants against the concentration gradient, i.e. a negative diffusivity (Deardorff and Willis, 1975). This feature may be

simulated by a higher-order closure turbulence model rather than the first-order closure assumptions of K,. (b) Vertical eddy diffusivity for neutral conditions. The formula given by Shir (1973) is adopted in the mode1for neutral conditions. This can be expressed as, K, = u,i, I= kzem4’lH

where H is the boundary layer height. (c) Vertical eddy diffusivity for stable conditions. The stable conditions, which typically occur at nighttime, exhibit weaker mixing than either the neutral or unstable cases. A profile given by Rao and Snodgrass (1979) for the steady-state nocturnal boundary layer is adopted in the model. K, =

ku,z 0.74 + 4.7 ZJL

ku,z 0.74+4.72/L

K, ;ii

Neutral

K, = ku,ze-4*/H

Unstable K, =

(2.20)

where b = 0.91, q = (z/f.)~- I”, and P = u,/vZ,I. lence closure model which was tested with stable boundary layer data from the Kansas and Minnesota experiments (Rao and Snodgrass, 1979).Further evaluation is ncccssary for application in an urban area under stable conditions when L is usually large or sometimes even negative due to the heat-island effect. The formulae for eddy difiusivities are summarized in Table 3. Figure 2 shows the vertical profiles of the eddy diffusivities in stable conditions and neutral conditions, and Fig. 3 shows the profile for unstable conditions. (d) Horizontal diffusion coeficients. In the horizontal direction, the transport due to turbulent diffu-

sion is much smaller than advection. The effective horizontal diffusion is thesum ofnumerical truncation

e-b”

152/L)*], 09 Z/Zi
e-br(

This equation was derived from a higher order turbu-

Table 3. Vertical eddy diffusjviti~ Stable

(2.19)

206

JIA-YEONG Ku er al.

for rtabkconditi~s and neutral conditions (if is the boundary height, u* is the frictional velocity and K, is the eddy diffusi~ty~.

Fig. 2. Vertical eddy diffusivity pro&

error ~nurnerj~l di~usion) and the true diffusion. The concentration disfritiution is found to be not very sensitive to the change of K, (&Rae et al., 1982). A value of SOm2 s - ’ is usedin the model (Reynolds er al,, 1973). 2.3.4. Dry depositionon the surface. The deposition rate is expressed as a function of the deposition velocity. A steady state balance is assumed between the turbulent transfer of pollutant down &concentration gradient and the net flux ofpollutant resulting from an exchange between the atmosphere and the earth’s surface. This is expressed by the lower boundary condition from Equation (2.2) as, - K*dCi/dZ = V&Ci(O) where C;(O) is the surface concentration and Uli is the deposition velocity of the species i which characterizes the interaction between the diffusing pollutant and the surface.

The specification of V, is difficult because of the complicated relationship of the deposition velocity to such factors as meteorological conditions, pollutant properties and surface characteristics. However, dry deposition rates should be proportjonal to the pollutant conctnfrafion. On this basis alone, one might expect dry deposition fluxes in urban areas to be somewhat greater than in the rural areas. The exchange of poJ&ants between the air and receptor surfaces is also IikeIy to he enhanced, ~~GWXZof the increased mechanical mixing caused by buildings and other structures in the urban areas. Sehmel (1980) reviewed deposition velocities based on field observations and laboratory experiments. The deposition velocity for SO2 ranged from 0.04 to 7.5 ems- ’ with values reported over St Louis among the highest. The deposition velocity for particulate sulfate is different from S02. SmaIl particles (for example, 0.1-i pm in diameter) are inefficiently deposited on

Numerical simuIation of air ~llution

207

in urban areas: modeldcve~opment

0.25

Fig. 3. Vertical eddy diffusivity profile for unstabk wnditions (Zi

is the. inversion height, we is the convective velocity scale and X, is the eddy dilkivity).

smooth surfaces (McMahon and Denison, 1979)due to smaller Brownian diffusion (compared to gases)and smaller inertial motion (compared to larger particles). The estimated base size distributions of sulfate particles in the St Louis area given by Pendergrass and Rao (1984)suggest that more than 70 % of the size of sulfate aerosol lies within the range 0.05-1 pm. In the model simulation, we assumed deposition velocities of Zcms- ’ and O.Zcms-’ for SO2 and sulfate, respectively. As a first approximation, the variation of deposition velocities as a function of meteorological conditions, including the stability class, is not considered in the model. 2.3.5. Chemical fransformation rare. In this model, chemical transformation of SO2 to sulfate is assumed to be linear and irreversible. The transformation rate depends on factors such as temperature, relative humidity and even the deposition velocity

(Alkezweenyand Powell, 1977).The average oxidation rate of SOz for the St Louis region varies from 2 to 16% h - ’ (Pendergrass and Rao, 1984).These studies were conducted both in summer and in fall seasons and the reported values regect daytime averages. A value of 5 Yeh- * is used in the model. For normal wind speeds( N 5ms-‘~it~~~that~er~tion~te is slower than the transport rate over a medium sized urban area. The time scaie of chemical transformation isT, = k;’ , where k. is the reaction rate constant. The advection time scale is T. = X/U where X is the diameter of the urban area, U is the mean wind speed. The chemical transformation will be important when U C X/T,. In St Louis X % 40 km which covers most of the measurement stations. If T, z 12 h, then U is about 1 ms" ‘. Thus, the time scales of chemical transformation and advection are comparable under low wind speed conditions.

JIA-YEONG Ku et al,

208

2.4 Methods of analysis 2.4.1. Grid specijcations. A three dimensional 30 x 40 x 7 grid system consists of 8400 grid elements (Fig. 4). The x, y, z axes, oriented along E-W, N-S, and vertical direction, are denoted by the indices j, k and I, respectively. The horizontal grid sizes are specified as follows:

Ax,=

l&Mm

6CjG25,

=2OOOm Ay, = 1OOOm =2OOOm

otherwise. 11
This “stretched” grid system covers an area of 40 x 60 km’. It encompasses smaller grid sub-elements compatible with the area source inventory grid size at

TUTM:725

the central urban area. The vertical grid sizes are specified as follows:

AZ, =

25 m, I = 1 5Om,2275m ! 50 m, 6 < I < 7 for H < 275 m.

The levels of the upper two grid points are determined by the hourly varying mixing depth, Zi. In this case, the top boundary H is always greater than or equal IO the height of the mixed layer. Allowance of grid points above Zi permits simulation of fumi~tion in the model. In the late morning hours when Z, rises rapidly, plumes trapped in the upper grid levels will be fumigated to the surface. 2.4.2. Numerical method. A time splitting method

XUTM-765

Fig. 4. The St Louis metropolitan area and horizontal numerical grid Bpacification

NumericA simulation of air pollution in urban areas: model development

209

as describedby Marchuk ff975) is used for integrating

the governing equation. The fourdimenrional partial differential equation given by (2.1)in (x, yVsl,t) is split into three tw~~rn~s~o~~ equations in (x, Q &, r)and kt). Each of these twodim~~ona~ ~~~tjons is integrated in succession over the same time step, and upon completion, the integration of the full fourd~~ona~ equation over one time step has been ap~rox~~t~. The three ~w~~i~n~~~~~~ equah3rrs

are:

A flux limitation (Rankoffand Hanzevack, 197.5)is imposed upon the system to prevent the co~ntratjon from becoming negative due to truncation error or phase error:

A similar formulation can be obtained in the y c~~~cnt. Xz/iit i 5/3X(&i) = ~~~~(~~~~j/~X) (22la) fb) Solution of the diffusive transport step+ The contribution to speciestransport from turbulent diffu(2.21b) BC;ISr+ S/Py(K,) = a/Jy(K,aCi/ay) sion depends on the coordinate direction. In the Howell piam_transport is dominated by adveetion Xji;it + Zj~iZfKj) = ~/~~~K~~~~~~~) + Q$i R and so, a simple, exphcit three-point finite difference i= I,2 (2*21c) can be employed for (T& and (T$.

solution can be obtained from the following sequence,

where n is the time level and Ar the time step, Each of the transport operators can he further @it into advtxtive and diffusive com~n~ts, ix 7, = U&V.&; G = (T,J.(T,b; and G = (K),, because of the assumption of zero vertical wind v&city. (a) Solution of the advective transport step. A second-order form~iat~onwritten in massconservative form is employed for the horizontal advection terms, (TX),and (T,),. These schemesare solved by an explicit method as follows:

where,

AXi-

r]/Axj,

(2.26)

An jrn~~~c~t ~n~te~i~ere~c~ ~~~~-~~c~olson method is used for (T,), to remove the time step limitation for the computational stability.

The numerical stability criteria in the x component knit the integration time step to he smaller of the advection limit%@At/Ax G $)+and the diffirsion limit (K,,Ar/Ax’ d i). The y component is done analogously. The diffusion term in the z component is solved by using an implicit finite difference technique to eliminate the s~~~~?y eon&aim. According to Richtmyer and Morton ff%?), the stability criteria for the implicit method are,

no restriction if 4 6 r d 1.

(2.29)

JIA-YEONG

210

This stability condition states that the implicit scheme is unconditionally stable when r 34. However, this cannot prevent amplification of truncation errors when the concentration fields are not smooth (large gradient) and/or the value of [At/(Az, + AZ, _ I)AzI] is large. The limitation is to make the diffusion flux vanish when the gradient of the concentration approaches zero. Thus, we have, (C :.;.‘I - C;,z.‘I+ 1MC;, II.I -q.

t., + I) 2 0.

This requirement leads to the condition K,,+iAt/(Az, + AZ, _ ,)Az, 6 1/[4( 1 -r)]. (2.30) Hence, the value r should be close to unity (0.98 is used in the numerical model simulations) when [At/(Az, + Azl _ ,)Az,] is large. The chemical transformation is limited by k.At C 1. 2.4.3. Poinr source/ormu/ulion. (a) Plume rise and effective release height. Plume rise calculation to determine the effective stack height (H,) is an important element in diffusion modeling. Except in strong winds, plume rise will increase the source height from 2 to 10 times of the actual release height H, for a typical elevated buoyant source. For most conditions, since the maximum ground level concentration is inversely proportional to the square of the effective stack height, this will reduce the maximum concentration by 3-100 times, or even more if the plume eflluent penetrates into an elevated stable layer (Briggs, 1984). Characterization of the plume rise above the stack top in terms of the exhaust gas properties (dominated by momentum or buoyancy) and the ambient atmospheric conditions is a complex problem. The reviews by Briggs (1975, 1984) indicate that no single formula can be applied to predict plume rise for the commonly encountered range of meteorological conditions. In convective or neutral conditions, the final plume rise is given

by either

of the two

following

equations

(Briggs, 1984). Ah = 30 (F/U)“:’

(2.3la)

Ah = 24 (F/U’)3” (H,+ 200 F/L”)““.

(2.3lb)

Ku et al.

rising into the mixed layer may partially or completely penetrate the elevated stable layer. Contributions to the ground level concentrations are from the fraction of the plume remaining within the mixed layer. This is more realistic than the “all or none” concept with respect to plume penetration of elevated inversions. The new penetration schemes used in this study follow the work of Briggs (1975, 1984), Weil and Brower (1984) and Rao (1985). The fraction “P” of the plume that penetrates the elevated stable layer is computed from Briggs’ (1984) criteria for bent over plumes. However, as a simple and conservative approach, Briggs assumes that the mixed-layer lapse rate is the same as that in upper stable layer, in order to account for the resistance of the inversion base to the plume penetration. The penetration criteria are based on the height of the upper and lower plume edges at final rise and the stable layer height above the stack Z: = Z, -H,, where Zi is the inversion height. Briggs (1975) assumed thal the heights of the upper and lower plume edges above the stack are about 1.5 and 0.5 times the plume rise, and the latter is given by: Ah’ = 2.6 (F/II Si)‘13

where si = (g/T,)d0,,/dz, and aO,,/dz is the potential temperature gradient in the stable layer. Then the penetration fraction P is defined as follows: P = 0,ifZ:/Ah’ > 1.5 P = 1,ifZi/Ah’ C 0.5 P = (1.5Ah’-Z;)/(l.SAh’-0.5Ah’),

Ah = 2.6 (F/Us)“’

(2.32)

and a9,/dz is the ambient where, s = (g/T,)dtl,/dr, potential temperature gradient. (b) Plume penetration into elevated stable layers. Quite frequently. the atmosphere consists of a well-mixed convective layer adjacent to the ground, capped by a stable layer at the top. A buoyant plume

(2.34)

if 0.5 < Z:,‘Ah’ < 1.5. The effective source height for plume material which is trapped in the mixed layer is given by Rao (1985) as: (2.35)

H; = H,++[Z:+Z:/(3-2P)].

The formulae for computing source strength and effective stack height are summarized as follows. If P = 1,the entire plume penetrates the stable layer and no ground level concentrations occur, i.e. H, 2 Zi and QS = 0.

Following a procedure described by Turner (1985) and Rao (1985), both equations are calculated and the lower value is used as the final plume rise. It is believed these equations represent a more realistic physical description of plume rise within the PBL. In stable condition, the final rise is given by:

(2.33)

(2.36a)

If P = 0, the plume remains within the mixed layer, i.e. H, = min (H,,H;) and QS= Q.

(2.36b)

If 0 < P < 1, the plume partially penetrates elevated stable layer, i.e. H,= H,+)(Z;+Z;/(3-2P))and

the

Q,= (1 -P)Q (2.36c)

where Q is the pollutant source strength at the stack exit and QS is the effective source strength of plume fraction remaining in the mixed layer. (c) Treatment of elevated point sources. One of the criticisms associated with grid models is related to the problem of accurately incorporating of the elcvatcd point sources into the grid system. One particular

Numerical simulation of air pollution in urban areas: model development

problem is related to the initial diffusion. However, once the plume size becomes comparable with the model grid scale, the behavior of the plume can be appropriately described by the atmospheric diffusion equation of the model. Therefore, the criterion for placing the plume material under control of the model must involve a comparison of the plume size with size of the grid cell in the model. The Lagrangian plume trajectory concept recommended by Shir and Shieh (1975) is employed in the model. The procedure keeps track of the location of the plume centerline and the plume width, and incorporates the plume into the grid system at the time when the plume width is comparable to the grid size. In other words, the effective point source location is calculated according to the wind speed and the diffusion speed. The method can be described as follows: r-=C+FT T = min (Ax*/K,,, Ay*/K,,, AZ*/&)

(2.37)

211

corporated in the lower boundary condition to treat the surface removal of the pollutants. A first-order chemical transformation process is also included for the treatment of reactive pollutants in the model. The performance of the numerical model formulated here is evaluated using the St Louis Regional Air Pollution Study (RAPS) data base in the following paper. Acknowledgemen&-Oneof the authors (J.-Y. Ku) would like to express his gratitude to the Atmospheric Turbulence and Dilfusion Laboratory, Oak Ridge, for providing the financial support, computer and library facilities and data. This work, which is part of his Ph.D. dissertation in the Department of Atmospheric Sciencesat the State University of New York at Albany (SUNYA). was accomplished under Interagency

Agreementsamong the U.S. Department of Energy, National Oceanic and Atmospheric Administration, the U.S. Environmental Protection Agency, New York State Department of Environmental Conservation and SUNYA Research Foundation.

(2.38)

where To,r’ and J are initial and new locations of point source, and wind vector at the initial point source location, respectively. The downwind distance at which plume material is mixed into the model grid cell depends on the rate of plume dispersion (i.e. on atmospheric stability). In the stable case, the plume will be dispersed slowly and could be carried out of the grid domain without being mixed into the model grid cells. The plume material is then mixed into a single cell, namely the cell that corresponds to the location of the plume centerline at the time the plume size equals the size of the grid cell. If the plume material were mixed into two or three grid cells over which a plume can actually extend, dispersion of the plume material would be overestimated, analogous to inducing an artificial dispersion for the plume.

3. SUMMARY

In this paper, the development of a numerical urban air pollution model which includes recent advances in treatment of pollutant transport and dispersion in the atmospheric boundary layer is described. Based on the mass conservation equation, the model incorporates advection, turbulent diffusion, chemical reaction, surface deposition and emissions of the pollutants, while reflecting the current understanding of planetary boundary layer dynamics to estimate and parameterize the atmospheric stability, eddy diffusivity and plume rise. The partial penetration of elevated stable layers by the plume is also considered in the model. A nondivergent wind field interpolation scheme which produces wind field with the structure inherent in the original wind measurements is included in the model in order to better represent the urban wind field characteristics and to maintain the physical constraint of mass conservation. The deposition velocity is in-

REFERENCE8 Alkezweeny A. J. and Powell D. C. (1977) Estimation of transformation rate of SO1 to SO, from atmospheric concentration data. A/mosphrric Environment I I, 179-182. BankolT S. G. and Hanzevack E. L. (1975) The adaptivefiltering transport model for prediction and control of pollutant concentration in an urban airshed. Atmospheric Enuironmenr9, 793-808. Blackadar A. K. (1962) The vertical diffusion of wind and turbulent exchange in a neutral atmosphere. J. Geophys. Res. 67, 3095-3 102. Briggs G. A. (1984) Plume rise and buoyancy effects. In Atmospheric Science and Power Production (edited by D. Randerson), Chapter 8. DOE/TIC-27601. Technical Information Center, Department of Energy, Oak Ridge, TN. Briggs G. A. (1975) Plume rise predictions. In Lectureson Air Pollurion and Environmental Impact Analysrs. pp. 59-l I I. Am. Meteor. Sot., Boston, MA. Businger J. A., Wyngaard J. C., lzumi Y. and Bradley E. F. (1971) Flux-profile relationships in the atmospheric surface layer. J. Ammos.Sci. 28, 181-189. Chang L.-F. W. (1978) Determination of surface flux of sensibleheat, latent heat, and momentum utilizing the bulk Richardson number. Papers in Meleorologirul Research 1. 16-24 (published in Taiwan). Clark T. L. and Eskridge R. E. (1977) Nondivergent wind analysis algorithm for the SI. Louis RAPS network. EPA600/4-77-049, EPA Research Triangle Park, NC. Crane G., Panofsky H. A. and Zeman 0. (1977) A model for dispersion from area sources in convective turbulence.

AtmosphericEnvironment 11,893-900. DeardorffJ. W. (1972)Numerical investigation of neutral and unstable boundary layer. J. Atmos. Sci. 29, 91-115. DeardorlTJ. W. and Willis G. E. (1975) A parameterization of dilTusion into the mixed layer. J. Appl. Meteor. 14, 1451-1458. Irwin J. S. (1979) A theoretical variation of the wind profile power-law exponent as a function of surface roughnessand stability. Atmospheric Enoironmenl 13, 191-194.

Kaimal J. S.,Wyngaard J. C, Haugen D.A., CotC0. R., ltumi Y.. Caughey S. J. and Readings C. J. (1976) Turbuknce structure in the convcclivc boundary layer. J. Armos. Sci. 33,2152-2169. Lamb R. G. ad Durran D. R. (1978) Eddy diffusivities derived from a numerical model of the convective boundary layer. I/ Nuouo Cirnenro lc, l-17 (published in Italy).

212

JIA-YEONC Ku er al.

Liu C. Y. and Goodin W. R. (1976) An iterative algorithm for objective wind field analysis. Mon. Wea. Rev. 104, 784792. Marchuk G. 1. (1975) Methods o/ Nmericol Mathemafics. Springer, New York. McMabon T. A. and Deaison P. J. (1979) Empirical atmospheric deposition parameters-a survey. Atmospheric Environment 13, 571-585. McRae G. J., Goodin W. R. and Seinfeld J. H. (1982) Development of a second-generation mathematical model for urban air pollution--I. Model formulation. Atmospheric Environment 16,679-6%. McRae G. J. and Seinfeld J. H. (1983) Development of a second-generation mathematical model for urban air pollution-H. Evaluation of model performance. Atmospheric Environment 17, 501-522. Pasquill F. (1962) Atmospheric Diflusion, p. 209. Van Nostrand, London. Pendergrass W. R. and Rao K. S. (1984) Evaluation of the Pollution Episodic Model using the RAPS data. NOAA Technical Memorandum ER< ARL-128, 47 pp. EPA600/3-84-087. U.S. EPA. Research Triannle Park. NC. Available from NTIS, P884-232537. Sprinifield, VA.~ _ Rao K. S. and SnodgrassH. F. (1979) Some parameterizations of Ihe nocturnal boundary layer. Boundary Layer Mereorology 17, 15-28. Rao K. S. (1985) User’s guide for PEM-2: Pollution Episodic Model (Version 2). NOAA Technical Memorandum ERLARL-137, 128 pp. To be published as a U.S. EPA report. Reynolds S. D., Roth P. M. and Seinfeld J. H. (1973) Mathematical modelingofpho~ochemical air pollution--l. Formulation of the model. Atmospheric Enuironmenf 7. 1033-1061. Richtmyer R. D. and Morton K. W. (1967) Diflerence

Methods/or Initial Value P&em% Inlerscience,New York. Schere K. L. (1981) Air quality model responseto objectively analyzed wind fields. Fi/lh Symposium on Turbulence, DiflLsion and Air Pollution. Atlanta, GA. Extended Abstracts, Amer. Meteor. Sot., Boston, pp. 42-43. Sehmel G. A. (1980) Particle and gas dry deposition: a review. Atmospheric Environment 14,983-1011. Shir C. C. (1973) A preliminary numerical study of atmospheric turbulent flows in the idealized planetary boundary iayer. J. Atnws. Sci. 30, 1327-1339. Shir C. C. and Shieh L. J. (19741 A generalized urban air pollution model and its applic&on-to Ihe study of SO* distributions in the St. Louis metropolitan area. J. Appl. Meteor. 13, 185-204. Shir C. C. and Shieh L. J. (1975) Development ofan Urban Air Quality Simulation Model wilh Compatible RAPS Data, Volume I. EPA-600/4-75-005a. U.S. EPA. Research Triangle Park, NC. Turner D. B. (1964) A diffusion model for an urban area. J. Appl.

Meteor.

3, 83-91.

Turner D. B. (1985) Proposed pragmatic methods for estimating plume rise and plume penetration through atmospheric layers. Atmospheric Environment 19, 12 I5- I2 18. Weil J. C. and Brower R. P. (1984) An update Gaussian plume model for tall stacks.J. Air Poll. Co&. Assoc. 34, 818-827. Willis G. E. and Deardorff J. W. (19781 A laboratorv studv of dispersion from an elevated’ so&e within a- modSled convective planetary boundary layer. Atmospheric Enuironmenl 12, 130>1312. Yamada T. (1977) A numerical experiment on pollutant dispersion in a horizontally-homogeneous atmospheric boundary layer. Atmospheric Environmen/ 11, 1015-1024. Yu T. W. (1977) A comparative study on parameterization of vertical turbulent exchangeprocesses.Mon. Wea. Rev. 105, 57-66.