# A Two-Dimensional Solution of the Advection-Diffusion Equation With Dry Deposition to the Ground

By Tirabassi, Tiziano Buske, Daniela; Moreira, Davidson M; Vilhena, Marco T

ABSTRACT A mathematical scheme is developed to simulate the vertical turbulent dispersion of air pollution that is absorbed or deposited to the ground. The scheme is an exact analytical solution of the atmospheric diffusion equation, without any restriction to the vertical profile of wind speed and eddy diffusivities, and taking into account the dry deposition by a boundary condition of a nonzero flux to the ground. The performances of the solution, with a proper parameterization of the vertical profiles of the wind and eddy diffusivities, were evaluated against the dataset from the Hanford (Washington) diffusion experiment, in which two tracers (one depositing and one nondepositing) were released simultaneously. In addition, the solution derived in this work is compared with four different models, with deposition at the ground, found in the literature.

(ProQuest: ... denotes formulae omitted.)

1. Introduction

The Eulerian approach for modeling the concentration of contaminants in a turbulent flow in the planetary boundary layer (PBL) is widely used in the field of air pollution studies. The simplicity of the K theory of turbulent diffusion has led to the frequent use of this theory as the mathematical basis for simulating air pollution dispersion. Dry deposition refers to the transfer of air pollution (gas and particles) to the ground, where it is removed. The various transfer mechanisms leading to dry deposition are complex and involve micrometeorological characteristics of the atmospheric surface layer (Seinfeld and Pandis 1998). The deposition flux is usually parameterized in terms of deposition velocity, which either is specified empirically or is estimated from appropriate theoretical relations. Moreover, the integral mass conservation equation is modified to account for the amount of material deposited on the surface.

The advection-diffusion equation can be written in finite- difference form, thus paving the way to a countless variety of numerical solutions. Using the gradient transport approach (K theory), dry deposition is included by specifying the deposition flux as the surface boundary condition. Therefore, numerical solutions to the advection-diffusion equation with variable eddy diffusivities are used to take into account the effects of dry deposition as well as gravitational settling for heavier particles (Arya 1999).

Analytical solutions of equations are of fundamental importance in understanding and describing physical phenomena (Pasquill and Smith 1983). Many operative models (using an analytical formula for the air pollution concentration) adopt empirical algorithms for describing dry deposition. The Gaussian plume equation was modified to include source depletion models (Chamberlain 1953; Overcamp 1976) and surface depletion models (Horst 1977, 1984). The solutions proposed by Smith (1962), Ermak (1977), and Rao (1981) also retained the framework of invariant wind speed and eddies with height (as the Gaussian approach). More recent analytical solutions of the advection-diffusion equation with dry deposition at the ground have utilized height-dependent wind speed and eddy diffusivities (Horst and Slinn 1984; Koch 1989; Chrysikopoulos et al. 1992; Lin and Hildemann 1997). However, these solutions are restricted to the specific case in which the source is located at the ground level and/ or with restrictions to the wind speed and eddy diffusivity vertical profiles.

We present an analytical solution of the advectiondiffusion equation obtained by means of the generalized integral Laplace transform technique (GILTT; Wortmann et al. 2005; Moreira et al. 2005), describing dry deposition with a boundary condition of nonzero flux to the ground and without any restriction to the above profiles and the source position. The GILTT has recently been applied for the simulation of pollutant dispersion in the atmosphere by solving analytically the advection-diffusion equation. The steps of this method are construction of an auxiliary Sturm-Liouville problem, expansion of the contaminant concentration in a series in terms of the obtained eigenfunctions, replacing of this equation in the original, and finally taking moments. The result is a set of ordinary differential equations that are then solved analytically by the Laplace transform technique. For more details, see the works of Wortmann et al. (2005) and Moreira et al. (2005).

We assume that concentration turbulent fluxes are proportional to the mean concentration gradient (first-order turbulence closure). The simplicity of the first-order turbulence closure (K closure) has led to its widespread use as the mathematical basis for simulating urban, photochemical pollution. However, K closure has its own limits. In contrast to molecular diffusion, turbulent diffusion is scale dependent. This means that the rate of diffusion of a cloud of material generally depends on the cloud dimensions and the intensity of turbulence. As the cloud grows, larger eddies are incorporated in the expansion process so that a progressively larger fraction of turbulent kinetic energy is available for the cloud expansion. However, eddies much larger than the cloud itself are relatively unimportant in its expansion. Thus, the gradient-transfer theory works well when the dimension of dispersed material is much larger than the size of turbulent eddies involved in the diffusion process (i.e., for ground-level emissions and for large travel times). In strict terms, one should introduce a diffusion coefficient function not only of atmospheric stability and emission height, but also of the travel time or distance from source. However, such time dependence makes it difficult to treat the diffusion equation in a fixed-coordinate system in which multiple sources have to be treated simultaneously. Otherwise, one should limit the application of the gradient theory to large travel times (Pasquill and Smith 1983). In addition, unlike molecular diffusion, turbulent diffusion is not a property of fluids but of the turbulence itself or of flows, and it may vary greatly from one flow to another and from one region to another of the same flow. The above relations are essentially based only on a qualitative analogy between molecular and turbulent diffusion. For the first-order closure to be realistic, the mean concentration field must have a much larger scale time than that of turbulent transport.

Despite these well-known limits, the K closure is widely used in several atmospheric conditions, because it describes the diffusive transport in an Eulerian framework, where almost all measurements are Eulerian in character (i.e., they are collected in a framework fixed with respect to the earth). It produces results that agree with experimental data as well as any more complex model, and it is not as computationally expensive as higher-order closures.

In this work, a step forward is taken toward solving the two- dimensional, steady-state advection-diffusiondeposition equation by the GILTT. The novelty depends on the construction of the auxiliary SturmLiouville problem. Indeed, for this kind of problem, the eigenvalues and eigenfunctions must be determined assuming boundary conditions of the third type, which encompass the contaminant deposition speed. Note that the mentioned works (Wortmann et al. 2005; Moreira et al. 2005) assume boundary conditions only of the second type. To validate the results obtained, numerical comparison is undertaken with available results in the literature.

2. The analytical solution

We present the construction of the two-dimensional analytical solution for the advection-diffusion-deposition equation to simulate pollutant dispersion in atmosphere with deposition to the ground, valid for any variable vertical eddy diffusivity coefficients and wind profile (without lateral dispersion). The concentration turbulent fluxes are often assumed to be proportional to the mean concentration gradient. This assumption, along with the equation of continuity, leads to the advection-diffusion equation. For a Cartesian coordinate system in which the x direction coincides with that of the average wind, the steady-state two-dimensional advection- diffusion equation with dry deposition to the ground is written as

... (1)

subject to the boundary conditions

... (1a)

... (1b)

and a continuous source condition

... (1c)

Here, C denotes the pollutant concentration, K1 is the turbulent eddy diffusivity coefficient assumed to be a function of the variable z, u is the mean wind oriented in the x direction and a function of the variable z, V^sub g^ is the deposition velocity, h is the height of PBL, Q is the emission rate, H^sub s^ is the height of the source, and delta is the Dirac delta function. The two- dimensional equation is equivalent to concentrations from an infinite line source or crosswind-integrated concentration from a point source.

To solve the problem by the GILTT, Eq. (1) is rewritten as

... (2)

where the first term on the right-hand side satisfies the following Sturm-Liouville problem: ... (3a)

... (3b)

The solution of the problem in Eq. (3) constitutes a well-known set of orthogonal eigenfunctions zeta^sub i^(z) = cos[lambda^sub i^(h - z)] whose eigenvalues fulfill the ensuing transcendental equation

... (3c)

where H^sub 1^ = Vg^sub l^[K^sub z^(z^sub 0^)]. [For more details, see the work of Ozisik (1980).] The eigenvalues are obtained from Eq. (3c) using the Newton-Raphson solving technique. The K^sub z^ is not evaluated at z = 0 but at z^sub 0^ (the roughness length) where V^sub g^ is defined also (Arya 1999). In fact, the lower boundary condition implies that the flux near the ground must be equal to the deposition rate, where also V^sub g^ is defined. Arya (1999), following van Ulden (1978), suggests z = z^sub 0^ as a good actual approximation of z - 0. Therefore, we followed the same approximation. Moreover, V^sub g^ is an empirical parameter and is obviously experimentally evaluated not at z = 0 but at some height from the ground, and it is considered to be constant until the ground is reached. Therefore, our approximation in the boundary condition follows the experimental approximation in the evaluation of V^sub g^.

It is now possible to apply the GILTT approach. For this purpose, the pollutant concentration is expanded in the series (Wortmann et al. 2005; Moreira et al. 2005)

... (4)

where zeta^sub i^(z) comes from the well-known solution of the Sturm-Liouville problem given in Eq. (3) and ... is given below.

By replacing the above equation in Eq. (1) and taking moments, the following is obtained (with j=1, 2, . . .):

... (5)

The above equation can be written in matrix fashion as

... (6)

where Y(x) is the column vector whose components are ..., the matrix F is defined as F = B-1E, and the entries of the matrices B and E are written by

... (7)

... (8)

Following the procedure of Wortmann et al. (2005) and Moreira et al. (2005), one obtains the following solution for the problem in Eq. (6):

Y(x) = XG(x)xi, (9)

where X is the eigenfunction matrix of F, G is the diagonal matrix whose entries have the form exp(-d^sub i^x), d^sub i^ are the eigenvalues of F, and xi is the vector given by xi = X^sup -1^Y(0). Applying the source condition in Eq. (1c) transformed by the GILTT technique, the unknown vector xi is determined by solving the resulting linear system. Therefore, the solution for the concentration given by Eq. (4) is now well determined because the vector ... is now known. For more details, see the work of Wortmann et al. (2005).

No approximations are made in the derivation of this solution, and therefore it is analytical except for the round-off error. The number of terms considered in the series summation is established in such a manner as to reach a prescribed accuracy, that is, a prescribed percent difference with the previous term of series summation (i.e., 0.5%). In this exercise, 30 terms were used. For better comprehension, Fig. 1 reports (as an example) the numerical convergence of the GILTT solution as a function of the number of eigenvalues.

To our knowledge, this is the first work that has solved the two- dimensional advection-diffusion equation by the GILTT method with a boundary condition of the third type, that is, nonzero flux at the ground. Observe that this solution is general in the sense that this solution for a boundary condition of the second type, that is, zero flux at the ground, is straightly attained by taking the limit when V^sub g^ [arrow right] 0 of the proposed solution. Here it is obvious that the eigenvalues and eigenfunctions changed, and, as before, they are tabulated in Ozisik (1980).

3. Experimental data and PBL parameterization

To show an example of the application of the obtained solution [Eq. (4)], the dataset of the Hanford diffusion experiment was used. This experiment was conducted in May-June 1983 on a semiarid region of southeastern Washington on generally flat terrain. The detailed description of the experiment was provided by Doran and Horst (1985). Data were obtained from six dual-tracer releases located at 100, 200, 800, 1600, and 3200 m from the source during moderately stable to near-neutral conditions. However, the deposition velocity was evaluated only for the last three distances. The release time was 30 min except in run 5, when it was 22 min. The terrain roughness was 3 cm.

Two tracers, one depositing and one nondepositing, were released simultaneously from a height of 2 m. Zinc sulfide (ZnS) was chosen for the depositing tracer, and sulfur hexafluoride (SF^sub 6^) was the nondepositing tracer. The lateral separation between the SF^sub 6^ and ZnS release points was less than 1 m. The near-surface release height and the atmospheric stability conditions were chosen to produce differences between the depositing and nondepositing tracer concentrations that could be measured easily. The data collected during the field tests were tabulated (as crosswind- integrated tracer concentration data) and were presented in Doran et al. (1984). The meteorological data and crosswind-integrated tracer concentration data, normalized by the release rate Q, are listed in Table 1. Note that in Table 1 C^sub d^ and C^sub nd^ are, respectively, the crosswind-integrated concentrations of ZnS and SF^sub 6^ normalized by the emission rate Q. For more details about the way that the effective deposition velocities and wind speed are calculated, see the work of Doran and Horst (1985). The PBL height h, not present in the Hanford dataset, was evaluated with the formula proposed by Zilitinkevich (1972):

... (10)

where f is the Coriolis parameter, u* is the vertical convective friction velocity scale, and L is the Monin-Obukhov length.

To use the above solution [Eq. (4)], it was necessary to select wind and eddy coefficient vertical profiles. The reliability of each model strongly depends on the way in which turbulent parameters are calculated and related to the current understanding of the PBL (Seinfeld and Pandis 1998).

The vertical eddy diffusivity used in this work is given in Mangia et al. (2002):

... (11)

where z is the height and Lambda = L(1 - z//z)^sup 5/4^.

The proposed K^sub z^ profile [discussed and derived in Mangia et al. (2002)] does not match the similarity profile for neutral conditions, as many other profiles proposed in the literature do not [e.g., the K^sub z^ vertical profiles proposed by Shir, Lamb et al., Businger and Arya, and Troen and Marth that are presented in Seinfeld and Pandis (1998)]. However, in the stable cases of the Hanford diffusion experiment, the profile of Eq. (11) follows the similarity profile in the surface layer, as we can see in Fig. 2.

The wind velocity profile was described by a power law expressed as follows (Panofsky and Dutton 1988):

... (12)

where u^sub z^ and u^sub 1^ are the mean wind velocity at the heights z and z^sub 1^ and n is an exponent that is related to the intensity of turbulence for rural terrain (Irwin 1979).

4. Numerical results

The model was evaluated with the ratio C^sub d^/C^sub nd^, where C^sub d^ and C^sub nd^ are the crosswind-integrated concentrations of ZnS and SF^sub 6^ measured at 1.5 m above the ground and normalized by the emission rate Q. A comparison of predicted and observed values C^sub d^/C^sub nd^ is shown in Fig. 3, with vertical eddy diffusivity given by Eq. (11) and the power profile of wind given by Eq. (12). In this respect, it is possible to note that the model reproduces fairly well the observed concentration.

Table 2 presents some performance measurements, obtained using the well-known statistical evaluation procedure described by Hanna (1989):

normalized mean square error (NMSE) = ...

factor of 2 (FA2) = fraction of data (%) within 0.5

correlation coefficient (COR) = ...

fractional bias (FB) = ... and

fractional standard deviation (FS) = ..., where subscripts 0 and p refer to observed and predicted quantities, respectively, sigma is the standard deviation, and an overbar indicates an average.

The statistical index FB indicates whether the predieted quantities underestimate or overestimate the observed ones. The statistical index NMSE represents the quadratic error of the predicted quantities in relation to the observed ones. The best results are indicated by values nearest 0 in NMSE, FB, and FS, and nearest 1 in COR and FA2. The statistical indices point out that a good agreement is obtained between experimental data and the GILTT model.

Doran and Horst (1985) presented four different models that evaluate the dry deposition at the ground with four different approaches: the source depletion approach of Chamberlain (1953), the corrected source depletion model of Horst (1980, 1983), the K model proposed by Ermak (1977) and Rao (1981), and the K corrected model of Rao (1981). For more information about the above models, see Doran and Horst (1985).

To make possible the comparison of our results with the ones in the models mentioned above, the statistical parameters used in the paper of Doran and Horst (1985) were adopted for the evaluation of their performances. These parameters are described in Fox (1981) and Willmott (1982):

mean bias ...

variance (...

mean absolute error ... and

index of agreement ...

where d^sub i^ is the difference between observed (C01) and predicted (C^sub pi^) values, ... the overbar indicates an average, 0

In Table 3, comparisons between the GILTT approach and the above models (Chamberlain 1953; Horst 1980,1983; Ermak 1977; Rao 1981) are reported, and it is possible to see the good performance of the solution presented here.

5. Final remarks

A new exact general solution of the two-dimensional steady state advection-diffusion equation has been presented that considers dry deposition to the ground and can be applied for describing the turbulent dispersion of many scalar quantities, such as air pollution, radioactive material, heat, and so on. To show the performances of the solution in actual scenarios, a parameterization of the PBL has been introduced, and the values predicted by the solutions have been compared with the Hanford diffusion experiment dataset. The analysis of the results shows a reasonably good agreement between the computed values and the experimental ones. The discrepancies with the experimental data depend not on the solution of the advection-diffusion equation but on the equation itself, which is only a model of the reality. Moreover, a source of discrepancies between the predicted and measured values lies in the PBL parameterization used (i.e., vertical wind and eddy diffusivity profiles). Last, the solution results were compared with those of four different models. Although models are sophisticated instruments that ultimately reflect the current state of knowledge on turbulent transport in the atmosphere, the results they provide are subject to a considerable margin of error. This is due to various factors, including in particular the uncertainty of the intrinsic variability of the atmosphere. Models, in fact, provide values expressed as an average, that is, a mean value obtained by the repeated performance of many experiments, whereas the measured concentrations are a single value of the sample to which the ensemble average provided by models refers. This is a general characteristic of the theory of atmospheric turbulence and is a consequence of the statistical approach used in attempting to parameterize the chaotic character of the measured data.

In light of the above considerations, an analytical solution is useful for evaluating the performances of sophisticated numerical dispersion models (which numerically solve the advection-diffusion equation), yielding results that can be compared not only with experimental data but, in an easy way, with the solution itself, to check numerical errors without the uncertainties presented above.

Acknowledgments. The authors thank the Conselho Nacional de Desenvolvimento Cientifico e Tecnologico (CNPq), the Coordenacao de Aperfeicoamento de Pessoal de Nivel Superior (CAPES), the Consiglio Nazionale delle Ricerche (CNR), and the project "Laboratorio LaRIA" for the partial financial support of this work.

REFERENCES

Arya, S., 1999: Air Pollution Meteorology and Dispersion. Oxford University Press, 310 pp.

Chamberlain, A. C., 1953: Aspects of travel and deposition of aerosol and vapour clouds. UKAEA Rep. AERE-HP/R-1261, Harwell, Berkshire, United Kingdom, 35 pp.

Chrysikopoulos, C. V., L. M. Hildemann, and P. V. Roberts, 1992: A three-dimensional steady-state atmospheric dispersion-deposition model for emissions from a ground-level area source. Atmos. Environ., 26A, 747-757.

Doran, J. C., and T. W. Horst, 1985: An evaluation of Gaussian plume-depletion models with dual-tracer field measurements. Atmos. Environ., 19, 939-951.

_____, O. B. Abbey, J. W. Buck, D. W. Glover, and T. W. Horst, 1984: Field validation of exposure assessment models. Data Environmental Science Research Laboratory, U.S. Environmental Protection Agency, Research Triangle Park, NC, EPA/600/384/092A, 177 pp.

Ermak, D. L., 1977: An analytical model for air pollutant transport and deposition from a point source. Atmos. Environ., 11, 231-237.

Fox, D. G., 1981: Judging air quality model performance. Bull. Amer. Meteor. Soc., 62, 599-609.

Hanna, S. R., 1989: Confidence limits for air quality model evaluations, as estimated by bootstrap and jackknife resampling methods. Atmos. Environ., 23, 1385-1398.

Horst, T. W., 1977: A surface depletion model for deposition from a Gaussian plume. Atmos. Environ., 11, 41-46.

_____, 1980: A review of Gaussian diffusion-deposition models. Atmospheric Sulfur Deposition, Environmental Impact and Health Effects: Proceedings of the Second Life Sciences Symposium, Potential Environmental and Health Consequences of Atmospheric Sulfur Deposition, Gatlinburg, Tennessee, October 14-18, 1979, D. S. Shriner, C. R. Richmond, and S. E. Lindberg, Eds., Ann Arbor Science, 275-283.

_____, 1983: A correction to the Gaussian source depletion model. Precipitation Scavenging, Dry Deposition and Resuspension, H. R. Pruppacher, R. G. Semonin, and W. G. N. Slinn, Eds., Elsevier, 1205- 1218.

_____, 1984: The modification of plume models to account for dry deposition. Bound-Layer Meteor., 30, 413-430.

_____, and W. G. N. Slinn, 1984: Estimates for pollution profiles above finite area-sources. Atmos. Environ., 18, 1339-1346.

Irwin, J. S., 1979: A theoretical variation of the wind profile power-low exponent as a function of surface roughness and stability. Atmos. Environ., 13, 191-194.

Koch, W., 1989: A solution of two-dimensional atmospheric diffusion equation with height-dependent diffusion coefficient including ground level absorption. Atmos. Environ., 23, 1729-1732.

Lin, J. S., and L. M. Hildemann, 1997: A generalized mathematical scheme to analytically solve the atmospheric diffusion equation with dry deposition. Atmos. Environ., 31, 59-71.

Mangia, C., D. M. Moreira, I. Schipa, G. A. Degrazia, T. Tirabassi, and U. Rizza, 2002: Evaluation of a new eddy diffusivity parameterisation from turbulent Eulerian spectra in different stability conditions. Atmos. Environ., 36, 67-76.

Moreira, D. M., M. T. Vilhena, T. Tirabassi, D. Buske, and R. M. Cotta, 2005: Near-source atmospheric pollutant dispersion using the new GILTT method. Atmos. Environ., 39, 6289-6294.

Overcamp, T. J., 1976: A general Gaussian diffusion-deposition model for elevated point sources. J Appl. Meteor., 15, 1167-1171.

Ozisik, M. N., 1980: Heat Conduction. John Wiley and Sons, 712 pp.

Panofsky, H. A., and J. A. Dutton, 1988: Atmospheric Turbulence. John Wiley and Sons, 397 pp.

Pasquill, F., and F. B. Smith, 1983: Atmospheric Diffusion. John Wiley and Sons, 437 pp.

Rao, K. S., 1981: Analytical solutions of a gradient-transfer model for plume deposition and sedimentation. NOAA Tech. Memo. ERL ARL-109, Air Resources Laboratories, 75 pp.

Seinfeld, J. H., and S. N. Pandis, 1998: Atmospheric Chemistry and Physics. John Wiley and Sons, 1326 pp.

Smith, F. B., 1962: The problem of deposition in atmospheric diffusion of paniculate matter. J. Atmos. Sci., 19, 429-434.

van Ulden, A. P., 1978: Simple estimates for vertical diffusion from sources near the ground. Atmos. Environ., 12, 2125-2129.

Willmott, C. J., 1982: Some comments on the evaluation of model performance. Bull. Amer. Meteor. Soc., 63, 1309-1313.

Wortmann, S., M. T. Vilhena, D. M. Moreira, and D. Buske, 2005: A new analytical approach to simulate the pollutant dispersion in the PBL. Atmos. Environ., 39, 2171-2178.

Zilitinkevich, S. S., 1972: On the determination of the height of the Ekman boundary layer. Bound-Layer Meteor., 3, 141-145.

TIZIANO TIRABASSI

Institute of Atmospheric Sciences and Climate of National Council of Researches, Bologna, Italy

DANIELA BUSKE

Institute of Atmospheric Sciences and Climate of National Council of Researches, Bologna, Italy, and Programa de Pos-Graduacao

em Engenharia Mecanica, Universidade Federal do Rio Grande do Sul, Porto Alegre, Brazil

DAVIDSON M. MOREIRA

Universidade Federal de Pelotas, UNIPAMPA, Bage, Brazil

MARCO T. VILHENA

Programa de Pos-Graduacao em Engenharia Mecanica, Universidade Federal do Rio Grande do SuI, Porto Alegre, Brazil

(Manuscript received 21 December 2006, in final form 15 January 2008)

Corresponding author address: Tiziano Tirabassi, Institute ISAC of CNR, Via Gobetti 101, 40129 Bologna, Italy.

E-mail: [email protected]

Copyright American Meteorological Society Aug 2008