A MODEL EMPLOYING INTEGRAL TRANSFORM METHOD TO SIMULATE POLLUTANT DISPERSION IN ATMOSPHERE

An updated version of the semi-analytical model for describing the steady-state concentration in the atmospheric boundary layer is presented here. Two inversion methods of the Laplace transform are tested: the Gaussian Quadrature scheme and the Fixed-Talbot method. The model takes into account settling velocity, removal (wet and dry deposition), and first order chemical reactions. The capability of the model to accurately predict the ground-level concentration is demonstrated qualitative and quantitatively. The results are in good agreement with experimental data.


INTRODUCTION
For many years, the atmospheric dispersion of pollutants has been rigorously studied in industrialized countries.The environmental problems caused by air pollution are complex and far-reaching, affecting many natural processes and strongly influencing the ecological balance.For this reason, it is important to improve our understanding of the dispersion process of pollutants in the atmosphere and the impact on the diverse ecosystems involved.Thus, the computational simulation of pollutant dispersion provides an important information source.Both our scientific understanding and technical developments have been greatly improved by the use of empirical, analytical, semi-analytical, and numerical models to predict the concentration variations in a plume.For this purpose, the advection-diffusion equation has been largely applied in operational atmospheric dispersion models.In principle, from this equation it is possible to obtain a theoretical model of dispersion from a source given appropriate boundary and initial conditions, plus knowledge of the mean wind velocity and concentration of turbulent fluxes.
Analytical solutions of advection-diffusion equations are of fundamental importance in describing and understanding dispersion phenomena (Pasquill and Smith, 1983), since all of the parameters are expressed in a mathematically closed form and therefore the influence of individual parameters on pollutant concentration can be easily examined.Also, the analytical solutions make it easy to obtain asymptotic behaviors of the solutions, which are usually difficult to obtain through numerical calculations.The analytical solutions can also be used to improve the modeling of pollutant dispersion (Ermak, 1977;Horst and Slinn, 1984;Koch, 1989;Chrysikopoulos et al., 1992;Lin and Hildemann, 1997;Arya, 1999;Wortmann et al., 2005;Moreira et al., 2005a;Szinvelski et al., 2006;Moreira et al., 2009;Cassol at al., 2009).In addition, the transport of pollutants in the atmosphere is characterized by turbulent diffusion, which has not been uniquely formulated in the sense that a single basic physical model capable of explaining all of the significant aspects of the transport process has not yet been proposed.
An alternative in the search of analytical solutions of differential equations is the use of the Laplace transform.Recognized by the scientific community, the Laplace transform is a powerful method for solving differential equations in engineering and science.The inverse of the Laplace transform is an important but computationally complex step in the application of the Laplace transform technique.The inverse Laplace transformation can be accomplished analytically, according to its definition, or by using Laplace transform tables.For a more complicated differential equation, however, it is difficult to analytically calculate the inverse Laplace transformation.Consequently, the numerical inverse Laplace transform algorithms are often used to calculate the numerical results.There are several numerical algorithms in the literature that can be used to perform the Laplace inversion (Davies and Martin, 1979;Narayanan and Beskos, 1982).Each individual method has its own application and is suitable for a particular type of function.
The focus of our paper is to present an updated semianalytical model for simulating pollutant dispersion in the atmospheric boundary layer (ABL) that considers two Laplace transform inversion methods and additional physical terms in the advection-diffusion equation.This semi-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 also in an easy way, with the solution itself to check numerical errors.The analytical feature and simplicity of the solution that will be shown in this work reinforces that the proposed method is a robust and promising method to simulate pollutant dispersion in atmosphere.
We examined two one-dimensional inversion routines: the Gauss Quadrature (GQ) algorithm and an updated version of the Fixed Talbot (FT) algorithm, which is based on deforming the contour in the Bromwich inversion integral.This approach is quite general in the sense it can be applied when the wind speed and eddy diffusivities are arbitrary continuous functions of height, but they are described by stepwise functions.For a better understanding of the mathematical features of the model, it is important to emphasize that in this approach no approximation is made in its derivation except for the stepwise approximation of the meteorological parameters and the Laplace numerical inversion.
To validate the results obtained, a numerical comparison is made with available datasets from the literature.To reach the objective, the paper is organized as follows: we start in section 2 by explaining the mathematical model with the Laplace inversion methods.In Section 3 we introduce the specific turbulent parameterizations used in the simulations.In Section 4 we briefly specify the experimental datasets.Then, in Section 5 we present numerical examples evaluating the performance of the two inversion algorithms and, finally, in Section 6 we draw conclusions.

THE MATHEMATICAL MODEL AND THE LAPLACE INVERSION METHODS
A typical problem with the advection-diffusion equation involves solving problems corresponding to continuous sources of pollution.More precisely, considering a Cartesian coordinate system in which the x direction coincides with that of the average wind, the steady two-dimensional advection-diffusion equation can be written as (Alam and Seinfeld, 1981): where C y is the crosswind integrated concentration, K z is the Cartesian component of eddy diffusivity in the vertical direction, u is the longitudinal wind speed, w is the vertical wind speed, w s is the setting velocity, α is the first order chemical reaction rate coefficient, and λ is the wet deposit.
The mathematical description of the dispersion problem represented by the Equation 1 is well posed when it is provided by boundary conditions.Indeed, it is assumed a source of constant emission rate Q: where ä (z -H s ) is the Dirac delta function and H s the source height.The pollutants are also subjected to the boundary conditions: (1) where h is the ABL height and V d is the dry deposition velocity.
In the following, K z , as well the wind speed u, depend only on the variable z and is assumed to be an averaged value.The stepwise approximation is applied in problem (1) by discretization of the height h into sub-layers in such a manner that inside each sub-layer, average values for K z and u are taken.It is important to note that this procedure transforms the domain of problem (1) into a multilayered-slab in the z direction.Furthermore, this approach is quite general in the sense it can be applied when these parameters are an arbitrary continuous function of the z variable.It is important to mention that u and K z are constant in each sub-layer, but the concentration still varies with x and z inside each layer.
It is now possible to recast problem (1) as a set of advective-diffusive problems with constant parameters, which for a generic sub-layer reads like: for n = 1:NL, where NL denotes the number of sub-layers and y n C the concentration at the n th sub-layer.Besides which, two boundary conditions are imposed at z = 0 and h (ABL height) given by Equation 2 together with the continuity conditions for the concentration and flux of concentration at the interfaces.Namely: must be considered, in order to uniquely determine the 2N arbitrary constants appearing in the solution of the set of problems (2).Now, applying the Laplace transform in Equation 2 results in: where , with source condition: and the boundary conditions: and which has the solution: where Applying the initial and boundary conditions, one obtains a linear system for the integration constants (A n and B n ).Then, the concentration is obtained by inverting numerically the transformed concentration C .At this point we need a method for the Laplace transform inversion.Because all methods have limitations, we emphasize the utilization of more than one algorithm to invert a transform.In the sequence we show the Laplace inversion methods.

Gaussian quadrature scheme
Firstly, the concentration is obtained by inverting numerically the transformed concentration by a Gaussian Quadrature (GQ) scheme: where a k and p k are the weights and roots of the Gaussian quadrature scheme tabulated in Stroud and Secrest (1966), M is the number of the quadrature points and, It is well known, that the results attained by the Gaussian quadrature scheme of order M are exact when the transformed function is a polynomial of degree (2N−1).It is possible, performing relative error calculations between the numerical results with k+1 and k points of quadrature, to control the error in the Gaussian quadrature scheme by properly choosing k, in order to attain a prescribed accuracy.

Fixed-Talbot scheme
Finally, a robust inversion method is applied, the Fixed Talbot (FT) algorithm (Talbot, 1979;Valkó and Abate, 2004;Abate and Valkó, 2004) and the concentration is obtained by: where and r is a parameter based on numerical experiments.The fixed Talbot (FT) method is based on the deformation of the contour of the Bromwich inversion integral and requires complex arithmetic.
Both GQ and FT methods have only one free parameter: M, which is the number of terms in the summation.Both algorithms provide increasing accuracy as M increases.Because all methods have limitations (see details see the works of Lee and Sheen (2004) and Moreira et al. (2005b)), we emphasize the utilization of more than one algorithm to invert a transform.
Since numerical Laplace inversion techniques are not exact, and often depend on the choice of a free parameter that is unknown a priori, it is advantageous to either use more than one inversion technique or perform experimentation and study the effect of the free parameter on the solution.In recent years, numerical transform inversion has become recognized as an important technique in operations research, notably for calculating probability distributions in stochastic and deterministic models.The significance of the numerical Laplace inversion is obvious from the wide range of applications.Well known in engineering, Laplace transformation methods are also used in order to solve differential and integral equations and to assist when other numerical methods are applied.
Laplace transforms are powerful tools used primarily to solve differential equations.The principal difficulty in using them is finding their inverses.Unless the transform is given in a table, an integration must be performed in the complex plane (Bromwich's integral) to find the inverse.Despite the power of complex analysis, this analytical technique often fails and Bromwich's integral must be integrated numerically.

TURBULENCE PARAMETERIZATIONS
In atmospheric diffusion modeling, the turbulent parameterization represents a fundamental aspect of the contaminant dispersion.The reliability of each model strongly depends on the way turbulent parameters are calculated and also on the current understanding of the ABL (Mangia et al., 2002).The literature reports many greatly varied formulae for the calculation of the vertical turbulent diffusion coefficient (Ulke, 2000).The cross-wind integrated concentrations from Prairie Grass experiment (Barad, 1958) and Hanford diffusion dataset (Doran and Horst, 1985) have been simulated using the formulae ( 8) and ( 10) and the following schemes for the estimation of eddy diffusivities and wind profile: Scheme I: employed in the convective conditions (Brost et al. 1988): where k is the von Karman constant (0.4), z is height above the ground, h is the thickness of the ABL and w * is the convective velocity.
Scheme II: employed in the stable conditions (Panofsky and Dutton, 1984): where u * is the friction velocity and the function Ö h is calculated with the formulae of Dyer: where L os the Monin-Obukhov lenght.
The equations used by the model to calculate mean wind are those predicted by similarity theory (Panofsky and Dutton, 1984): where u is the scale velocity relative to mechanical turbulence, z 0 is the roughness length and m Ø the stability function expressed in Businger et al. (1971) relations: ( ) (stable conditions) and (convective conditions) where

EXPERIMENTAL DATA AND NUMERICAL RESULTS
The performance of the Eulerian semi-analytical model has been evaluated against experimental ground-level concentrations provided by the Prairie Grass (Barad, 1958) and Hanford (Doran et al. 1984) diffusion experiments.

Comparison with Prairie Grass data det -unstable case
The Prairie Grass experiment was conducted in O'Neill, Nebraska in 1956.The pollutant (SO 2 ) was emitted without buoyancy at a height of 0.5 m and it was measured by samplers at a height of 1.5 m in five downwind distances (50, 100, 200, 400, 800 m).The Prairie Grass site was flat with a roughness length of 0.6 cm.The results for twenty convective (-h/L > 10) experiments are presented.For more details see the work of Nieuwstadt (1980).
Table 1 presents some performance measurements, obtained using the well-known statistical evaluation procedure described by Hanna (1989).The statistical index FB indicates whether the predicted 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.
In general, the concentrations from scheme I are close to those observed within a factor of 2 of 74% and correlation of 94% for Laplace inversion of FT with M = 20.The scheme GQ with M = 20 predict the concentrations with a factor of 2 of 71% and correlation of 92%.The confidence interval for FB and NMSE calculated (not shown) partly intersect each other, showing that there is not certainty that the results of the model are different for any scheme of inversion.
A comparison of predicted and observed values is shown in Figure 1, with vertical eddy diffusivity given by scheme I and the profile of wind given by Equation 15.In this respect, it is possible to note that the model reproduces fairly well the observed concentration using FT inversion method.
In addition, we deduce from Figure 1 that the model is able to more accurately reproduce the concentrations closer to the source (higher concentrations), where turbulence is based on superficial scales.The discrepancies increase at farther distances from the source (lower concentrations), where the plume penetrates the free convection layer, rapidly decreasing the concentration (Nieuwstadt, 1980).

Comparison with Hanford data set -stable case
The Hanford experiment was conducted in May-June 1983 on a semi-arid region of southeastern Washington on generally flat terrain.A detailed description of the experiment is provided by Doran and Horst (1985).Data were obtained from six dual-tracer releases located at 100, 200, 800, 1600, and    16) 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 non-depositing, were released simultaneously from a height of 2 m.Zinc sulfide (ZnS) was chosen for the depositing tracer, and sulfur hexafluoride (SF6) was the non-depositing tracer.The lateral separation between the SF6 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 non-depositing tracer concentrations that could be measured easily.The data collected during the field tests were tabulated (as crosswindintegrated tracer concentration data) and were presented in Doran et al. (1984).For more details about the way that the effective deposition velocities and wind speed are calculated, see the work of Doran and Horst (1985).Table 2 presents model performance results with Hanford experimental data for ZnS.
The concentrations from scheme II are close to those observed within a factor of 2 of 100% and correlation of 92% for Laplace inversion of FT with M = 20 (number of terms in summation) considering deposition.The schemes FT with M = 20 (non-depositing) and GQ with M = 20 (depositing) predict the concentrations within a factor of 2 of 22% and a correlation of 83 and 85%, respectively.The results for these two last situations are not favorable.
A comparison of predicted and observed values is shown in Figure 2, with vertical eddy diffusivity given by scheme II and the profile of wind given by Equation 15.It is possible to note that the model reproduces very well the observed concentration when considering the boundary condition (5c) with the FT method.
The effect of the boundary condition (5c) is illustrated in Figure 2. Figure 2 compares the concentrations both with and without deposition and the corresponding depletion of material from the plume.A pollutant with zero deposition velocity (white circles) is totally reflected back into the atmosphere.By contrast, a pollutant with the highest deposition velocity (black squares) undergoes substantial depositional losses, thereby reducing the airborne concentrations.Increasing the dry deposition reduces the extent of human exposure.At low emission heights, pollutant removal by dry deposition should explicitly be taken into account.Especially under stable atmospheric conditions, the plume is restricted to a narrow surface layer and the interaction of the plume with the surface can cause considerable mass removal (Koch, 1989).
It is possible to observe a large difference in the results of Tables 1 and 2 between the methods FT and GQ.In Table 1 there is good agreement between the simulation results using FT and GQ.However, in Table 2, there is a large discrepancy between the methods.The following Figures may explain this , and ).We observed in Figure 3 that for high sources the discontinuity existing due to the presence of the source using GQ method deteriorates the solution.This fact does not occur with FT method.However, the concentration at ground level is coincident for both methods (in fact, the more important concentration is at ground level).In Figure 4 we see that the presence of the source close to the ground for the GQ method (unlike the FT method) did not produce good results.To more accurately assess these points shown in Figures 3 and 4, we show below simulations of ground-level concentration (Cy/Q (10 -3 s.m -2 )) using the data from experiment 1 of Hanford (hh = 325m, L = 165m, u = 3.6m/s, u * u * = 0.4m/s).
Figure 5 shows a good convergence with increasing parameter M using the FT method for a low source of 3m.However, in Figure 6 we observe that even with the increased number of terms M in the GQ method the convergence does not occur, indicating that for a low source (close to the ground) the     method don't have good results (including negative values of concentration for distances greater than 1000m for M > 2).In Figure 6, distances less than 1000m don't show negative values of concentration for M > 8, which makes the results for M = 20 in the GQ method very similar to FT in Table 1.The methods show very similar results for ground-level concentration with M = 20 for distances less than 1000m.
Figures 7 and 8 show the results for a source of 100 m height.It is possible to observe in Figure 7 a good convergence with increasing M in the GQ method.In Figure 8, we present results from using the FT method.We observed oscillations for distances close to the source, but with the increase of M the solution has good convergence with M 20.It should be noted that these simulations are for concentrations at ground level.The source in Figures 5 and 6 was located near the ground.
Such methods are susceptible to the Gibbs phenomenon near discontinuities and this is a problem particularly when working with distributions defined on a finite range if the function does not decline smoothly to zero at the end-points.Gibbs phenomena are small oscillations around the edges of sharp gradients.Depending on the system these oscillations will, in the worst case, eventually deteriorate or destroy the solution.
We observe in the simulations differences in the results when using FT and GQ methods.This analysis was done for the first time and shows that for low sources the most suitable method is FT.This methodology was used for atmospheric problems, but can be extended to other areas of knowledge because the method is computationally simple and robust.

CONCLUSIONS
A solution of the two-dimensional steady state advectiondiffusion equation has been presented that considers dry and wet deposition, settling velocity, first order chemical reactions, and can be applied for describing the turbulent dispersion of many scalar quantities, such as air pollution, radioactive material, and heat.This solution is general for a boundary condition of the second type, that is, zero flux at the ground is straightly attained by taking the limit when V d = 0 in the boundary condition of the proposed problem.To show the performances of the solution in actual scenarios, two Laplace transform inversions have been used and the values predicted by the solution have been compared with the Prairie Grass and Hanford diffusion experiments dataset.The analysis of the results shows very good agreement between the computed values and the experimental ones when using a Laplace inversion with the FT method.This study shows that among these two methods (GQ and FT), the FT transform inversion technique is the most powerful but also the most computationally expensive.The method can be used in many engineering applications and is easy to implement and leads to accurate results for many problems including diffusion-dominated and other functions.The GQ method fails to predict the concentration for low sources at long distances.Such methods are susceptible to the Gibbs phenomenon near discontinuities and this is a problem particularly when working with distributions defined on a finite range if the function does not decline smoothly to zero at the end-points.Gibbs phenomena are small oscillations around the edges of sharp gradients.Depending on the system these oscillations will in worst case eventually deteriorate or destroy the solution.
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 ABL parameterization used (i.e., vertical wind and eddy diffusivity profiles).
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.The aptness of the method shows that it is capable of solving realistic problems.Moreover, the analytical feature and simplicity of the solution reinforces that the proposed method is a robust and promising method to simulate pollutant dispersion in atmosphere.

Figure 1 -
Figure 1 -Scatter diagram of observed and predicted data.Data between the middle diagonal line indicates perfect agreement.Dotted lines indicate a factor of two.

Figure 2 -
Figure 2 -Scatter diagram of observed and predicted data.Data between the middle diagonal line indicates perfect agreement.Dotted lines indicate a factor of two.

Figure 5 -
Figure5-Ground-level concentration as a function of source distance using FT method in Laplace inversion for a source height of 3m.

Figure 6 -
Figure6-Ground-level concentration as a function of source distance using GQ method in Laplace inversion for a source height of 3m.

Figure 7 -
Figure7-Ground-level concentration as a function of source distance using GQ method in Laplace inversion for a source height of 100m.

Figure 8 -
Figure8-Ground-level concentration as a function of source distance using FT method in Laplace inversion for a source height of 100m.

Table 1 -
Statistical evaluation of model results for Prairie Grass experiment.

Table 2 -
Statistical evaluation of model results for Hanford experiment.