Variable K-theory model for the dispersion of air pollutants in low-wind conditions in the surface-based inversion

ARTICLE IN PRESS

Atmospheric Environment 41 (2007) 6951–6963 www.elsevier.com/locate/atmosenv

Variable K-theory model for the dispersion of air pollutants in low-wind conditions in the surface-based inversion Maithili Sharan, Manish Modani Centre for Atmospheric Sciences, Indian Institute of Technology, Delhi, Hauz Khas, New Delhi 110016, India Received 14 May 2006; received in revised form 31 October 2006; accepted 13 November 2006

Abstract A variable K-model has been proposed for the dispersion in low winds in the surface-based inversion by expressing the eddy diffusivities as a linear function of downwind distance from the source. The resulting partial differential equation with variable coefficients along with the physically relevant boundary conditions is solved analytically. For the accuracy of the so-obtained solution, an analysis of the convergence and error estimation has been carried out. It is shown that the series converges absolutely. An upper bound for the error based on the partial sum of the series is estimated and it is described that the error tends to zero as the number of terms in the expansion are sufficiently large. The solution has been used to simulate the field tracer data sets collected from Hanford and IIT diffusion experiments in stable and unstable conditions, respectively. It predicts 41% cases in stable and 35% cases in unstable conditions within a factor of two to observations. r 2006 Elsevier Ltd. All rights reserved. Keywords: Mathematical model; Variable eddy diffusivities; Eigen-function expansion; Low wind dispersion; Inversion/mixing layer

1. Introduction The commonly used mathematical models to estimate pollutant concentrations, are known to work in all meteorological conditions except low winds. In these conditions, the gaps in the dispersion models appear to exist at two levels: (i) at the formulation level that involves certain approximations/assumptions and (ii) at the level of applications, for example, nonavailability of suitable dispersion parameters. In most of the classical models, the down-wind diffusion is neglected in Corresponding author. Tel.: +91 11 26591312; fax: +91 11 26591386. E-mail address: [email protected] (M. Sharan).

1352-2310/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.atmosenv.2006.11.018

comparison to the advection and this may incur an error upto 25% (Sharan et al., 1996a; Sharan and Modani, 2005). Thus, it has been emphasized that the inclusion of longitudinal diffusion in the models would be more appropriate for the treatment of dispersion in low winds (Arya, 1995). Over the years, attempts have been made to develop air quality models for the dispersion of air pollutants in low-wind conditions accounting for the longitudinal diffusion (Sirakov and Djolov, 1979; Arya, 1995; Sharan et al., 1996a, b) etc. However, these models account for the unrestricted diffusion in the vertical direction. This is in contrast to the fact that the atmospheric boundary layer (ABL) is often capped by an inversion which tends to reflect back the material hitting the inversion base (Arya, 1999).

ARTICLE IN PRESS 6952

M. Sharan, M. Modani / Atmospheric Environment 41 (2007) 6951–6963

Recently, Sharan and Modani (2005) described a mathematical model accounting for the diffusion in all-three directions with the advection along the mean wind and restricted vertical diffusion. The eddy diffusion coefficients are assumed constant to obtain the analytical solution of the resulting partial differential equation. However, for practical application, the so-obtained solution was expressed in terms of downwind distance using the traditional relationship between eddy diffusion coefficients and dispersion parameters. The approach to derive the solution assuming eddy diffusion coefficients constant and later on, introducing the dependence in terms of downwind distance is mathematically inconsistent (Llewelyn, 1983). Moreira et al. (2005) also presented a steady state mathematical model by solving three-dimensional advection–diffusion equation using Laplace transform and considering the planetary boundary layer as a multilayer system. In each sublayer, the averaged value of eddy diffusion coefficients was assumed. The model was validated for peak concentrations only in both stable and unstable conditions. However, it is known that the longitudinal diffusion term does not contribute to the peak concentrations occurring along the plume centerline. Researchers are continuously attempting to evolve a general parameterization scheme for eddy diffusivities to develop more realistic and accurate models for the dispersion of pollutants in the atmosphere especially in low-wind conditions. Arya (1995) emphasized that the eddy diffusivities should be parameterized as a function of downwind distance for the treatment of near-source dispersion. Later, Sharan et al. (1996b) proposed a model by parameterizing eddy diffusion coefficients in terms of downwind distance using Taylor’s (1921) theory of plume dispersion. However, this model accounts for the unrestricted diffusion in the vertical direction. Though the model proposed by Sharan and Modani (2005) deals with the dispersion in a finite layer, the eddy diffusivities are assumed constant in deriving the solution of the resulting partial differential equation. In addition, limitations are associated with its applications in terms of dispersion parameters in low-wind conditions. In order to overcome some of these, a more realistic model is proposed for the treatment of dispersion in low-wind conditions in a surface-based inversion. The model is validated with Hanford and IIT diffusion data sets in stable and unstable conditions, respectively.

2. Mathematical description The commonly used advection–diffusion equation in K-theory-based air-quality models for steady-state concentration c in the absence of removal (physical or chemical) mechanisms and orienting x-axis in the direction of mean wind U is given as       qc q qc q qc q qc U ¼ Kx Ky Kz þ þ , qx qx qx qy qy qz qz (1) where Kx, Ky and Kz are the eddy diffusion coefficients along x, y and z directions, respectively. Arya (1995) has suggested simple parameterization of eddy diffusivities based on Taylor’s (1921) statistical theory and similarity relations for use in a variable K-theory (gradient-transport) model. He has argued that for small travel times (i.e., small distances from the source) diffusivities can be expressed as a linear function of downwind distance; the constant of proportionality represents the turbulence parameters. This parameterization also accounts for atmospheric stability through the turbulence parameters, which essentially represent turbulent intensities. Therefore, here, the eddy diffusion coefficients for nonzero wind speed are parameterized as K x ¼ aUx;

K y ¼ bUx;

K z ¼ gUx ,

(2)

where a, b and g represent the turbulence intensities. Assuming U constant, Eq. (1) reduces to ax

q2 c qc q2 c q2 c þ bx þ ða  1Þ þ gx ¼ 0. qx2 qx qy2 qz2

(3)

Eq. (3) is subject to the following boundary conditions: (i) A continuous point source with strength q is assumed to be located at the point (0, 0, H), i.e. Uc ¼ qdðyÞdðz  HÞ at

x ¼ 0,

(4)

in which d(.) is the Dirac’s delta function and H the effective stack height. (ii) Concentration of the pollutant tends to zero far away from the source, i.e.   c ! 0 as x;  y ! 1. (5)

ARTICLE IN PRESS M. Sharan, M. Modani / Atmospheric Environment 41 (2007) 6951–6963

(iii) The pollutant is not absorbed by the ground surface and therefore, there is no diffusive flux at the surface, i.e. K z

qc ¼0 qz

at

z ¼ 0.

(6)

(iv) The pollutant dispersion remains confined to a layer capped by an inversion lid at the top, which serves as an impermeable upper boundary for the pollutant (Dobbins, 1979; Beyrich, 1997). Thus, the boundary condition at the top (h) of the inversion/mixed layer can be prescribed as

6953

2.1. Special cases Case I: When source is located at the ground: The expression for the concentration is obtained from the relation (8) by taking H-0 and it comes out to be   2 m   q 1 x 1 pffiffiffiffiffiffi cðx; y; zÞ ¼ G m þ Uh GðmÞr bp ar2 2 1 1

 pffiffiffiðmþ12Þ X rn mþ2 þ4 r g 2 n¼1

 pffiffiffi    K mþ12 rn r g cos rn z . ð9Þ

An analytical solution of Eq. (3) for nonzero a with the boundary conditions (4)–(7) is obtained using methods of eigen-function expansion and Fourier transforms (Appendix A):   2 m   q 1 x 1 pffiffiffiffiffiffi cðx; y; zÞ ¼ G m þ Uh GðmÞr bp ar2 2 1 1

 pffiffiffiðmþ12Þ X rn mþ2  pffiffiffi þ2 r g K mþ12 rn r g 2 n¼1

      cos rn ðz þ H Þ þ cos rn ðz  H Þ ,

Case II: When there is no inversion (i.e. h-N): This case essentially represents the dispersion of a pollutant in the absence of inversion lid (Hanna et al., 1982). The limit h-N in Eq. (8) has been obtained using the property of Riemann integral to convert the summation series into an integral (Sharan and Modani, 2005). The resulting expression can be written as q pffiffiffiffiffi cðx; y; zÞ ¼ 2 2px U bg "  ðmþ1Þ a y2 ð z  H Þ 2 þ  1þ 2 x b g   ðmþ1Þ # a y2 ð z þ H Þ 2 þ þ 1þ 2 . x b g

ð8Þ

ð10Þ

in which m ¼ 1/2a, Km is the modified Bessel function of second kind of order m, rn ¼ np/h and r2 ¼ ðx2 =aþ y2 =bÞ. The ratio p ¼ (x2/ar2) indicates the closeness of the receptor to the plume centerline. The value of p ¼ 1 implies that the receptor lies on the plume centerline whereas po1 represents the extent of the closeness of the receptor to the plume centerline. In other words, p closer to one represents the region in the vicinity of the plume centerline. Notice that c decreases as p decreases implying that the concentration decreases away from the plume centerline. The solution (8) represents the distribution of pollutant’s concentration emitted from an elevated point source in presence of the inversion/mixed layer and accounting for the diffusion in all-three coordinate directions. It may be interpreted as the contribution from a source term located at z ¼ H and its successive images by reflection from two parallel boundaries at z ¼ 0 and z ¼ h (Hanna et al. 1982; Seinfeld, 1986). Notice that the concentration distribution represented by (8) is not Gaussian.

Solution (10) is the same as obtained by Sharan and Yadav (1998) for the dispersion of a pollutant in the absence of surface-based inversion.

K z

qc ¼0 qz

at

z ¼ h.

(7)

Case III: Slender plume approximation (i.e. a-0): This approximation gives the concentration close to the plume centerline (Sharan et al., 1996b) and amounts to neglect the downwind diffusion term in Eq. (1). The resulting expression is obtained by solving Eq. (3) for a ¼ 0 following the approach as discussed for a6¼0 in Appendix A, i.e.   q y2 pffiffiffiffiffiffiffiffi exp  cðx; y; zÞ ¼ 2bx2 Uhx 2pb "   1 X x2 g 2 r  1þ exp  2 n n¼1 #   cosðrn ððz þ HÞÞ þ cosðrn ðz  HÞÞ . ð11Þ

ARTICLE IN PRESS M. Sharan, M. Modani / Atmospheric Environment 41 (2007) 6951–6963

6954

Replacement of bx2 by sy2 and gx2 by sz2 in the above equation yields the Gaussian plume solution ! q y2 cðx; y; zÞ ¼ pffiffiffiffiffiffi exp  2 2sy 2pUsy h "   1 X s2  1þ exp r2n z 2 n¼0 #   cosðrn ðz þ HÞÞ þ cosðrn ðz  HÞÞ , ð12Þ where sy and sz are the dispersion parameters in y and z directions, respectively.

3. Results and discussion For the accuracy of the obtained solution (8), an analysis of the convergence and error estimation has been carried out (Appendix B). The series converges absolutely. The upper bound for the error tends to zero as the number of terms in the expansion are sufficiently large for all x and z. This is also confirmed computationally as Fig. 1 shows that the concentration becomes almost constant as number of terms in the expansion increases. In order to examine the performance of the present model, concentrations predicted by the model have been compared with those observed during the Hanford and IIT diffusion experiments in stable and unstable conditions, respectively.

3.1. Stable conditions The Hanford diffusion experiment was conducted in May–June 1983 at Hanford on a semi-arid region of south eastern Washington (461340 N, 1191360 W) on generally flat terrain. The detailed description of the experiment was provided by Doran and Horst (1985). Totally 6 test runs were conducted (Table 1). Samplers were placed at angular distances of 81, 41, 41, 21 and 31 on the concentric circles (center at release point) of radii 100, 200, 800, 1600 and 3200 m, respectively. The release height of SF6 was 2 m and average release rate was around 0.3 g s1. Meteorological inputs have been provided by the measurements taken from a tower located at 100 m to the north of the release point (Doran et al., 1984). The values of parameters such as atmospheric stability, U, root mean square of horizontal (sy) and vertical (sf) wind directions, u* and h are given in Table 1. The terrain was considered as an urban terrain (Doran and Horst, 1985). Concentrations have been computed from the present model (Eq. (8)). The parameters a, b and g are identified as squares of turbulence intensities using Taylor’s statistical theory of diffusion (Arya, 1995): a¼

s 2 u

U

;

b¼

s 2 v ; U

g¼

s 2 w

U

120000 Concentration

100000 80000 60000 40000 20000 0 50

(13)

where su, sv and sw are the standard deviations of velocity components in the longitudinal, lateral and vertical directions, respectively, and in general, these are available from the turbulence measurements. In the absence of these measurements, the su and sv are estimated from the following formulations

140000

0

,

100

150

200

250

n Fig. 1. Concentration computed from Eq. (8) with the number of terms for x ¼ 100 and z ¼ 0.5.

ARTICLE IN PRESS M. Sharan, M. Modani / Atmospheric Environment 41 (2007) 6951–6963

6955

Table 1 The meteorological data observed at 2 m height during the Hanford diffusion grid experiment Run no.

Date (DD/MM/YY)

P-G stability

U (m s1)

sy (m)

sf (m)

u* (m s1)

h (m)

1 2 3 4 5 6

18-05-83 26-05-83 05-06-83 12-06-83 24-06-83 27-06-83

E E E E E E

3.63 1.42 2.02 1.50 1.41 1.54

13.0 12.8 13.1 13.7 15.8 14.9

4.1 3.5 3.7 3.8 4.6 4.7

0.40 0.26 0.27 0.20 0.26 0.30

325 135 182 104 157 185

(Hanna et al., 1982): z

su ¼ 2u 1  , (14) h z

sv ¼ 1:3u 1  . (15) h sw is calculated from the relation (Arya, 1999) sw ¼ tanðsf Þ, (16) U where sf is in radians. For small values of sf, tan(sf)Esf and accordingly Eq. (16) becomes sw ¼ sf . (17) U This approximation is, in fact, used by Leahey et al. (1988) and Weber (1998). As the values of sf in both these studies are less than 61 whereas these values are less than 51 (Table 1) in the present study, the formulation (17) is applicable in present case as well. Diffusion data have been analyzed for overall concentrations by restricting computations upto the arcs of 800 m radii as our interest is to treat nearsource dispersion only. As an illustration, the observed and computed concentrations obtained from present model for test run #3 are shown in Fig. 2. In this test run, present model overpredicts on 100 m arc, predicts comparable concentrations on 200 m arc and gives underpredicting trend on 800 m arc in comparison to those observed. The scatter diagram (Fig. 3(a)) between the ratios of predicted to observed concentrations and observed concentrations shows that about 40% cases are predicted within a factor of two to those observed. Specifically, present model predicts 41%, 68%, 74%, 79% and 80% cases within factors of 2, 4, 6, 8 and 10 to those of observations, respectively. For the sake of comparison, the concentrations obtained from the Constant-K model developed by

Sharan and Modani (2005) are also shown in Figs. 2 and 3(a). The dispersion parameters s’s required in Constant-K model are estimated from the empirical formulations proposed by Briggs (1973). The constant-K model predicts 38%, 56%, 62%, 68% and 68% cases within factors of 2, 4, 6, 8 and 10 to the observations, respectively. The statistical measures such as fractional variance (FS), geometric mean bias (MG) and geometric mean variance (VG) (Hanna et al., 1991) have been computed to analyze the performance of the present model quantitatively. Ideal values for FS, MG and VG are 0, 1 and 1, respectively. Magnitude of MG less than one indicates overprediction whereas greater than one represents underprediction of the concentrations. The values of FS, MG and VG obtained from the present model are 0.64, 0.47 and 95.92, respectively. It shows that present model overpredicts. However, the values of FS and MG for Constant-K model are 0.09 and 13.07, respectively. It shows that Constant-K model relatively underpredicts. In case of Constant-K model, VG could not be estimated because at all the points beyond the range of 73sy the magnitude of the computed concentrations is low. The FS, MG and VG for concentrations predicted in the range of 73sy from Constant-K model are 0.02, 0.95 and 4.05, respectively. The unpaired analysis, i.e. Quantile–Quantile (Q–Q) plot (Olesen, 1995; Venkataram, 1999) (Fig. 3(b)) shows that the concentrations predicted for the observations larger than 250 ppt are within a factor of two to observations from the present model. Below this level, the present model overpredicts. However, in case of Constant-K model, concentrations predicted for the observations larger than 750 ppt are within a factor of two to observations and below this, it severely underpredicts. This may be due to the fact that Constant-K model with Briggs formulation is

ARTICLE IN PRESS M. Sharan, M. Modani / Atmospheric Environment 41 (2007) 6951–6963

6956

70000

Concentration (ppt)

60000 50000 40000 30000 20000 10000 0 18000 Concentration (ppt)

16000 14000 12000 10000 8000 6000 4000 2000 0 1600 Concentration (ppt)

1400 1200 1000 800 600 400 200 0 80

90

100 110 120 Angular distance (degrees)

130

Fig. 2. Observed and predicted concentrations for run # 3 on (a) 100, (b) 200 and (c) 800 m. Observed ( ( ) and IGP (JJJJ).

known to predict concentrations at the points lying in the range of 73sy. The overpredicting trend of present model may be due to the uncertainties associated with the parameters such as h, u*, sf, etc. The results from the present model are also compared with those computed from a steady-state model obtained by integrating the Gaussian puff (IGP) formula in finite vertical domain with respect to time between 0 and N, by considering the ~ and sz ¼ g~ t. The parameterization sx ¼ a~ t; sy ¼ bt

140

), present (

), constant-K

expression for c is given as follows: cðx; y; zÞ ¼

 rffiffiffi  2 2 q0 u2 p 1 u x pffiffi exp exp  2 ~ 2 r~ 2~a 2~a4 r~ 2ph~ab   X 1    ux cos rn ðz þ H Þ  erfc  pffiffiffi 2 pffiffi þ 2a~ r~ n¼1 Z   1 1 þ cos rn ðz  H Þ t2 0   g~ 2 t2 1 ux dt ,  exp  2 r~ þ 2  rn ð18Þ 2t 2 a~ t

ARTICLE IN PRESS M. Sharan, M. Modani / Atmospheric Environment 41 (2007) 6951–6963

6957

Predicted/Observed Concentrations

100 Present Constant-K

10

1

0.1

0.01

0.001 10

100

1000 Observed Concentrations

10000

100000

100000 Present Constant-K

Predicted Concentrations

10000 1000 100 10 1 0.1 0.01 1

10

100 1000 Observed Concentrations

10000

100000

Fig. 3. (a) Ratio of predicted and observed concentrations with observed concentrations in Hanford diffusion experiment. Middle line is the one-to-one line and the outer lines are the lines with factor two. (b) Q–Q plot between the observed and predicted concentrations in ppt in Hanford diffusion experiment. Middle line is the one-to-one line and the outer lines are the lines with factor two.

where erfc is the complementary error function and ! x2 y2 þ r~ ¼ . a~ 2 b~ 2 A model on the similar arguments in the absence of inversion has been discussed by Cirillio and Poli (1992). The parameters a~ ; b~ and g~ are obtained from the following formulations (Cirillio and Poli, 1992): a~ ¼ su ;

b~ ¼ sv ;

g~ ¼ sw .

(19)

The parameters su, sv and sw have been estimated from Eqs. (14), (15) and (17), respectively. Fig. 2

reveals that both the present and IGP models give rise to almost similar results. The value of FS, MG and VG are obtained as 0.67, 0.41 and 151.9, respectively, which are also close to that obtained from the present model. This is also expected as both the models have similar basic structure and use the same type of parameterizations (Sharan and Yadav, 1998). 3.2. Unstable conditions In unstable conditions, the data obtained from IIT-SF6 diffusion experiment, conducted at IIT Delhi (281520 N, 771180 E) in low-wind conditions

ARTICLE IN PRESS M. Sharan, M. Modani / Atmospheric Environment 41 (2007) 6951–6963

6958

during February 1991, has been used. Details of the experiment including atmospheric stability are given in Singh et al. (1991) and Sharan et al. (1996a). The release rate of SF6 tracer varied between 30 and 50 ml min1. Wind and temperature measurements were obtained at four levels (2, 4, 15 and 30 m) from 30-m micrometeorological tower. In the absence of onsite sounding data, the mixing heights were obtained from sodar measurements taken by National Physical Laboratory at a nearby location in Delhi (Sharan et al., 1996b). The values of wind speed, atmospheric stability, mixing height and u* are given in Table 2. Layer-averaged mean wind has been used while simulating the concentration distribution at the sampling height z ¼ 0.5 m. Table 2 reveals that the atmospheric conditions during the experiment were of weakly unstable type. The sufficiently large values of u* and U (Table 2) exclude the possibility of the free-convection case. Therefore, in the absence of direct measurements, the dispersion parameters a, b and g are computed

Table 2 Meteorological conditions during the IIT diffusion experiment Run no.

Time

U (m s1)

h (m)

Atmospheric stability

u* (m s1)

1 2 6 7 8 11 12 13

1200–1230 1530–1600 1000–1030 1245–1315 1645–1715 1000–1030 1215–1245 1530–1600

1.56 0.74 1.34 1.54 0.89 1.07 1.55 1.08

1570 1240 1070 1240 943 1070 1325 1070

A–B B B B B A–B B B

0.34 0.21 0.34 0.37 0.25 0.25 0.35 0.21

in terms of friction velocity (Sharan et al., 2002), i.e. u 2  a ¼ 4:28 ; U

b ¼ 4:0

u 2 

U

and

g ¼ 1:69

u 2 

U

.

(20) The peak concentrations computed from Eq. (8) on 50 and 100 m arcs are tabulated along with the observations in Table 3. It shows that the peak values computed from the present model are within a factor of 2 to observations in most of the runs. In run 7, the observed concentration at 50 m is suspect because it is lower that at 100 m arc. For the sack of comparison, table also includes the results computed from the Constant-K model and Briggs formulations (Briggs, 1973), respectively. The results indicate that the simulated peaks from Constant-K model are better than those obtained from the present model. The NMSE for the peak concentrations from the present and Constant-K models are 0.38 and 0.18, respectively. It is noticed that NMSE for the peak concentrations, computed from a steady-state model proposed by Moreira et al. (2005) is 0.35, which is close to that obtained from the present model. The scatter diagram for the overall concentrations (both on and off centerline) from the present and Constant-K models are given in Fig. 4(a). The figure reveals that the results from the present and Constant-K models are comparable. The cases predicted in a factor of two to observations are 35% from both the models and the corresponding values of NMSE are 1.41 and 1.07, respectively. The Q–Q plots (Fig. 4(b)) shows that the concentrations predicted for the observations larger than 5 ppt are within a factor of two to observations

Table 3 Observed and predicted peak values of tracer concentrations (ppt) at 50 and 100 m downwind distance Run no.

50 m arc Observed

1 2 6 7 8 11 12 13

832 1068 1101 248 1282 616 759 1060

100 m arc Predicted

Observed

Present

Constant-K

444 668 444 349 877 486 574 783

1025 1012 405 345 1015 656 752 1447

345 460 176 288 345 162 222 215

Predicted Present

Constant-K

137 167 111 87 219 122 144 196

318 252 102 87 256 165 189 364

ARTICLE IN PRESS M. Sharan, M. Modani / Atmospheric Environment 41 (2007) 6951–6963

6959

Predicted/Obsserved Cocnentrations

100 Present Constant-K

10

1

0.1

0.01 1

10

100 Observed Concentrations (ppt)

1000

10000

Predicted Concentrations (ppt)

10000 1000

Present Constant-K

100 10 1 0.1 0.01 1

10

100 Observed Concentrations (ppt)

1000

10000

Fig. 4. (a) Ratio of predicted and observed concentrations with observed concentrations in IIT diffusion experiment. Middle line is the one-to-one line and the outer lines are the lines with factor two. (b) Q–Q plot between the observed and predicted concentrations in IIT diffusion experiment. Middle line is the one-to-one line and the outer lines are the lines with factor two.

from the present model. Below this, present model underpredicts. The underpredicting trend of the present model in majority of the cases (Figs. 4(a) and (b)) may be due to the uncertainties associated with the measurement of u*, used in the parameterization of su, sv and sw, respectively. We wish to point out that slight decrease (upto 0.05 m s1) in the values of u* improves the performance of the present model and reduces the NMSE from 1.41 to 0.83. 4. Conclusions A variable-K mathematical model accounting for the diffusion in the all-three coordinate directions

and advection along the mean wind in a finite inversion layer has been proposed for the treatment of pollutant’s dispersion under low-wind conditions. The eddy diffusivities have been parameterized in terms of downwind distance for near-source dispersion. An analytical solution for the resulting advection–diffusion equation with the physically relevant boundary conditions has been obtained using the method of eigen-function expansion and Fourier transform. The particular cases corresponding to various physical situations from the soobtained analytical solution have been deduced. For the accuracy of the solution, an analysis for the convergence and error estimation has been carried out.

ARTICLE IN PRESS M. Sharan, M. Modani / Atmospheric Environment 41 (2007) 6951–6963

6960

The model is validated with the data sets obtained from Hanford diffusion experiment in stable conditions and IIT diffusion experiments in unstable conditions, respectively. The performance of the present model has been compared with the Constant-K model proposed by Sharan and Modani (2005) in both stable and unstable conditions, respectively. In addition, the concentrations obtained from the present model in stable condition have also been compared with the model obtained from the integration of Gaussian puff formula in finite layer with respect to time from 0 to N. Therefore, further validation of the present model is desirable with the availability of diffusion data in elevated release under low-wind conditions. The upstream diffusion may be expected near the source under low-wind conditions (Arya, 1995). However, the solution obtained in the present study has a practical limitation that it does not give the concentration field in the region upstream to the source. The parameterization of eddy diffusivities used in this study is valid for nonzero wind speed. In case of calm wind (i.e. zero wind speed) and unrestricted diffusion of pollutants in the vertical directions, the parameterization suggested by Arya (1995) can be used. However, a general parameterization of eddy diffusivities for the treatment of near-source dispersion is desired for the formulation of a more realistic dispersion model valid for all wind speeds.

A.1. Estimation of unknown coefficients Differentiate Eq. (A1) with respect to x, y and z and substituting them in Eq. (3) to obtain  axA0xx þ ða  1ÞA0x þ bzA0yy 1  X þ axAnxx þ ða  1ÞAnx n¼1

þbzAnyy  gxr2n An cosðrn zÞ ¼ 0,

ðA2Þ

where the suffixes x and y denote the partial derivatives with respect to x and y, respectively. Using orthogonal properties of eigen-functions (i.e. multiplying Eq. (A2) by cos rjz, j ¼ 0; 1; . . . and integrating it from 0 to h with respect to z), we get a set of partial differential equations: axA0xx þ ða  1ÞA0x þ bzA0yy ¼ 0,

(A3)

axAnxx þ ða  1ÞAnx þ bzAnyy  gxr2n An ¼ 0, n ¼ 1; 2; . . .

ðA4Þ

The boundary conditions (Eqs. (4) and (5)) become A0 ð0; yÞ ¼

qdðyÞ Uh

An ð0; yÞ ¼

2q dðyÞ cosðrn HÞ Uh

at x ¼ 0,

(A5)

n ¼ 1; 2; 3; . . .

at x ¼ 0

ðA6Þ

and An ðx; yÞ ! 0;

Acknowledgments We gratefully acknowledge valuable comments and suggestions from the anonymous reviewers. Authors also wish to thank Dr. Anil Kumar Yadav for his comments on the manuscript.

Appendix A In view of the boundary conditions (Eqs. (6) and (7)), the solution of Eq. (3) can be written as cðx; y; zÞ ¼ A0 ðx; yÞ þ

1 X

An ðx; yÞ cosðrn zÞ,

(A1)

n¼1

where {A0, A1, A2, y} is the set of unknown coefficients. Here rn ¼ np/h, n ¼ 0, 1, 2, y are the eigen-values and cos (rnz) are the corresponding eigen-functions.

x; jyj ! 1;

n ¼ 0; 1; 2; . . .

(A7)

Notice that Eq. (A3) with boundary condition (A5) may be deduced from Eq. (A4) and the corresponding boundary condition (A6) by taking n ¼ 0 and q in place of 2q. Thus, we will present briefly the solution of Eq. (A4) with the boundary conditions (A6) and (A7) only. Use the Fourier transform with respect to y, Z þ1 A~ n ðx; lÞ ¼ An ðx; yÞe2pily dy, (A8) 1

in Eq. (A4) to obtain an ordinary differential equation   axA~ nxx þ ða  1ÞA~ nx  4p2 l2 b þ gr2n xA~ n ¼ 0. (A9) The solution of the ordinary differential equation (A9) is a modified Bessel function of order m ¼ 1/2a provided a6¼0. Its complete solution can

ARTICLE IN PRESS M. Sharan, M. Modani / Atmospheric Environment 41 (2007) 6951–6963

be given by 8 9 1 = <  2 b g A~ n ðx; lÞ ¼ xm 4C n I m x 4p2 l2 þ r2n : ; a a 8 93 1 = <  2 b g 5, þDn K m x 4p2 l2 þ r2n : ; a a 2

ðA10Þ in which Cn and Dn are unknown coefficients. Since xm Im (x)-N and xm Km(x)-0 as x-N, for a bounded solution Cn ¼ 0. The unknown coefficient Dn is estimated by applying the transformed condition U A~ n ð0; lÞ ¼ ð2q=hÞ cos rn H corresponding to Eq. (A6). Using the approximation of Bessel function Km(X) ¼ (G(m)/2)(X/2)m for small argument (Gradshteyn and Ryzhik, 1980) we get Dn ¼

22m q sinðmpÞGð1  mÞ p Uh  m g 2 2 2 2b  4p l þ rn cosðrn HÞ. a a

ðA11Þ

22m q sinðmpÞGð1  mÞ A~ n ðx; lÞ ¼ p Uh  m g 2 2 m 2 2b x 4p l þ rn a a 8 9  1 = < 2 b g cosðrn HÞ.  K m x 4p2 l2 þ r2n : ; a a ðA12Þ Inverting (A12) with respect to l using d 6.726.4 (Gradshteyn and Ryzhik, 1980), we find  m pffiffiffimþ12 g 4 q b An ðx; yÞ ¼ pffiffiffi sinðmpÞGð1  mÞ r a p Uh r mþ12  pffiffiffi  n x2m K mþ12 rn r g cosðrn HÞ, 2 ðA13Þ

2

r ¼



 x2 y2 þ . a b

Once Ai’s are obtained, solution c in Eq. (A1) is given by   2 m   q 1 x 1 p ffiffiffiffiffi ffi cðx; y; zÞ ¼ G mþ Uh GðmÞr bp ar2 2 1 1  pffiffiffi mþ1 X rn mþ2  pffiffiffi þ 2 r g ð 2Þ K mþ12 rn r g 2 n¼1

      cos rn ðz þ H Þ þ cos rn ðz  H Þ . ðA14Þ Appendix B. Convergence and error analysis The expression (Eq. (8)) for concentration c(x, y, z) is rewritten as 1 mþ1 X  pffiffiffi rn 2 K mþ12 rn r g cðx; y; zÞ ¼ f 0 ðx; yÞ þ f ðx; yÞ 2 n¼1       cos rn ðz þ H Þ þ cos rn ðz  H Þ , ðB1Þ

Substituting Cn and Dn in (A10), we obtain

where

6961

where

   m q G m þ 12 x2 pffiffiffiffiffiffi , f 0 ðx; yÞ ¼ Uh GðmÞr bp ar2  pffiffiffiðmþ12Þ  m q 2 r g x2 pffiffiffiffiffiffi f ðx; yÞ ¼ , Uh GðmÞr bp ar2  2  np x y2 and r2 ¼ rn ¼ þ . h a b

B.1. Convergence The absolute value of c is written as  1  mþ1    X  pffiffiffi  rn 2      jcjp f 0 ðx; yÞ þ 2 f ðx; yÞ K mþ12 rn r g .  2 n¼1 P1

(B2)

Using the result, n¼1 jan j converges, P if the series a converges as well. For this, then the series 1 n n¼1 we make use of D’ Alembert test (Krasnov et al., 1990):  pffiffiffi    mþ12 K mþ12 rnþ1 r g anþ1  r nþ1  ¼ Lt  pffiffiffi . Lt  n!1 an  n!1 rn K mþ12 rn r g (B3) For large n, using the pffiffiffiffi ffi approximation of Bessel p expðX Þ for large argument function K m ðX Þ ¼ 2X

ARTICLE IN PRESS M. Sharan, M. Modani / Atmospheric Environment 41 (2007) 6951–6963

6962

(Gradshteyn and Ryzhik, 1980) to get    m pffi pffi anþ1   ¼ Lt 1 þ 1 ephr g ¼ ephr g . Lt   n!1 an n!1 n p pffi h40; r40 and g40 implies ehr g o1.

the remaining terms in the asymptotic expansion vanish for the large argument. (B4) (B5)

This implies that the infinite series (B2) of absolute terms converge; indicating the infinite series (B1) also converges. B.2. Error analysis We make use of the concept of partial sum of the series for error estimation. Define the partial sum of the series (B1) as m mþ1 X  pffiffiffi rn 2 Sm ¼ f 0 ðx; yÞ þ f ðx; yÞ K mþ12 rn r g 2 n¼1       cos rn ðz þ H Þ þ cos rn ðz  H Þ and accordingly the error     jE m j ¼ S mþp  S m p2 f ðx; yÞ  mþp X  r mþ12  pffiffiffi   n . K 1 rn r g  mþ2  2 

ðB6Þ

n¼mþ1

Let N be a fixed number sufficiently large such that m+pXN for all p. Since the series (B1) is convergent, (m+1)th term is smaller than the mth term. This implies amþp oaN

for all p.

Based on Eq. (B7), Eq. (B6) simplifies to      rN mþ12  pffiffiffi     jE m jp2 f ðx; yÞ p K mþ12 rN r g . 2

(B7)

(B8)

Use the approximation of Bessel function for large argument as N become large and substitute f (x, y) to get pffiffiffi 2 m  m pffi Npr g gx 2q p Np  h pffiffiffi jE m jo e . (B9) Uh GðmÞr b ar2 2h For each m, the error Em is bounded. As mO(1) and the remaining variables on the right hand side of Eq. (B9) are finite, the upper bound tends to 0 as N-N because the exponential term decreases faster than Nm. We wish to point out that we have used an approximation of Bessel function for large argument in Eq. (B8) by retaining the first term in the asymptotic analysis of the Km (Watson, 1995). However, our conclusion will not be affected as

References Arya, S.P., 1995. Modeling and parameterization of near-source diffusion in weak winds. Journal of Applied Meteorology 34, 1112–1122. Arya, S.P., 1999. Air Pollution Meteorology and Dispersion. Oxford University Press, New York, 320pp. Beyrich, F., 1997. Mixing height estimation from Sodar data—a critical discussion. Atmospheric Environment 31, 3941–3953. Briggs, G.A., 1973. Diffusion estimation for small emissions. US National Oceanic and Atmospheric Administration E.R.L. Report ATDL-106. Cirillio, M.C., Poli, A.A., 1992. An inter-comparison of semiempirical diffusion models under low wind speed, stable conditions. Atmospheric Environment 26A, 765–774. Dobbins, R.A., 1979. Atmospheric Motion and Air Pollution. Wiley, New York, 343pp. Doran, J.C., Horst, T.W., 1985. An evaluation of Gaussian plume-depletion models with dual-tracer field measurements. Atmospheric Environment 19, 939–951. Doran, J.C., Abbey, O.B., Buck, J.W., Glover, D.W., Horst, T.W., 1984. Field validation of exposure assessment models. EPA/600/384/092A, vol. 1, Data Environmental Science Research Lab, NC, 177pp. Gradshteyn, I.S., Ryzhik, I.M., 1980. Table of Integrals, Series and Products. Academic Press, New York, 1160pp. Hanna, S.R., Briggs, G.A., Hosker Jr., R.P., 1982. Handbook on atmospheric diffusion. US Department of Energy Report COE/TIC-11223, Washington, DC, 102pp. Hanna, S.R., Strimaitis, D.G., Chang, J.C., 1991. Evaluation of commonly-used hazardous gas dispersion models. In: Hazard Response Modeling Uncertainty. A Quantitative Method, vol. II, AFESC Contract No. FO8635-89-C-0136, H.Q. AFESC/RDVS, Tyndall AFB, FL 32403, 1991. Krasnov, M., Kiselev, A., Makarenko, G., Shikin, E., 1990. Mathematical Analysis for Engineers. Vol. I, Mir Publishers, Moscow. Leahey, D.M., Hansen, M.C., Schroeder, M.B., 1988. An analysis of wind fluctuation statistics collected under stable atmospheric conditions at three sites in Alberta, Canada. Journal of Applied Meteorology 27, 774–777. Llewelyn, R.P., 1983. An analytical model for the transport, dispersion and elimination of air pollutants emitted from a point source. Atmospheric Environment 17, 249–256. Moreira, D.M., Tirabassi, T., Carvalho, J.C., 2005. Plume dispersion simulation in low wind conditions in stable and convective boundary layers. Atmospheric Environment 39, 3643–3650. Olesen, H.R., 1995. The model validation exercise at Mol: overview of results. International Journal of Environment and Pollution 5, 761–784. Seinfeld, J.H., 1986. Atmospheric Chemistry and Physics of Air Pollution. Wiley, New York, 738pp. Sharan, M., Modani, M., 2005. An analytical study for the dispersion of pollutants in a finite layer under low wind conditions. Pure and Applied Geophysics 162, 1861–1892.

ARTICLE IN PRESS M. Sharan, M. Modani / Atmospheric Environment 41 (2007) 6951–6963 Sharan, M., Yadav, A.K., 1998. Simulation of diffusion experiments under light wind, stable conditions by a variable K-theory model. Atmospheric Environment 32, 3481–3492. Sharan, M., Singh, M.P., Yadav, A.K., Agarwal, P., Nigam, S., 1996a. A mathematical model for dispersion of air pollutants in low wind conditions. Atmospheric Environment 30, 1209–1220. Sharan, M., Singh, M.P., Yadav, A.K., 1996b. A mathematical model for the atmospheric dispersion in low winds with eddy diffusivities as linear functions of downwind distance. Atmospheric Environment 30, 1137–1145. Sharan, M., Yadav, A.K., Modani, M., 2002. Simulation of short-range diffusion experiment in low-wind convective conditions. Atmospheric Environment 36, 1901–1906. Singh, M.P., Agarwal, P., Nigam, S., Gulati, A., 1991. Tracer experiments—a report. Tech. Report, CAS, IIT Delhi.

6963

Sirakov, D.E., Djolov, G.D., 1979. Atmospheric diffusion of admixtures in calm conditions. Geophysique 32, 891–892. Taylor, G.I., 1921. Diffusion by continuous movements. Proceedings of the London Mathematical Society, Series 2 XX, 196–212. Venkataram, A., 1999. Applying a framework for evaluating the performance of air quality models. In: Proceedings of the Sixth International Conference on Harmonisation Within Atmospheric Dispersion Modelling for Regulatory Applications, Rouen, France, 11–14 October 1999. Watson, G.N., 1995. A Treatise on the Theory of Bessel Functions, second ed. Cambridge Mathematical Library, 812pp. Weber, R.O., 1998. Estimators for the standard deviations of lateral, longitudinal and vertical wind components. Atmospheric Environment 32, 3546–3639.