An Analytical Dispersion Model for Sources in the Atmospheric Surface Layer with Dry Deposition to the Ground Surface

This study presents a dispersion model based on an analytical solution of two-dimensional advection-diffusion equation to simulate the crosswind integrated concentration in the bounded atmospheric boundary layer. The dispersion model considers the power law profiles of wind speed and vertical eddy diffusion coefficient and deposition velocity at the ground surface. A closed form analytical solution of the advection-diffusion equation for these parameterizations of wind and eddy diffusivity was derived by considering the deposition through lower boundary condition at the ground surface. A sensitivity analysis of the dispersion model for ground level concentrations with varying deposition velocities at the ground surface was performed. Hanford diffusion experiment for a depositing tracer was used to evaluate the dispersion model in stable condition. Model was evaluated by accounting both deposition and non-deposition at the ground surface. The simulated concentrations from the proposed dispersion model with dry deposition are in good agreement with those observed and the model predicts ~72% cases within a factor of two. It was shown that the consideration of no-deposition condition in the model for a depositing tracer gives rise to a larger over-prediction to the observations in comparison to with deposition.


INTRODUCTION
Dry deposition is an imperative removal process for some air pollutants in the atmosphere to study and alleviate the air pollution problems in the atmosphere. The influence of dry deposition process on the turbulent dispersion of a passive contaminant in the atmospheric boundary layer (ABL), through the gradient transport theory based advectiondiffusion equation (Seinfeld, 1986), is usually described by taking into account the deposition flux as a surface boundary condition. An analytical dispersion model is always of great importance to fast and efficiently analyzing the effect of an individual input parameter on the simulated results.
Gaussian models that considered the normal distribution of concentrations in vertical and lateral directions are widely accepted and have been used to study the dispersion process in the atmosphere and air quality analysis for regulatory purposes. However, Gaussian models have some intricacy to take account the effect of particle settling, deposition and the spatial variability of the meteorological variables on the dispersion (Smith, 2008). To account the removal process in Gaussian models, an analytical solution of the advectiondiffusion equation with constant wind speed and eddy diffusivity is generally used (Overcamp, 1976;Horst, 1977Horst, , 1984. However, the assumption of invariant profiles of wind speed and eddy diffusivity in the solution of advectiondiffusion equation is not physically realistic and it limits the applicability of these modified Gaussian dispersion models. This assumption of constant wind speed and eddy diffusivity was further released over the years by considering them as the functions of space coordinates in derivation of the solution of the advection-diffusion equation (Chrysikopoulos et al., 1992;Lin and Hildemann, 1997;Essa et al., 2007;Vilhena et al., 2008;Moreira et al., 2009;Sharan, 2010, 2012). Chrysikopoulos et al. (1992) and Lin and Hildemann (1997) derived the exact solutions of the advection-diffusion equation with dry deposition on the ground surface and for power law profiles of the vertical eddy diffusivity and wind speed in the unbounded atmosphere (infinite mixing/inversion layer) for the ground level area and point sources, respectively. However, an assumption of infinite unbounded ABL in derivation of these solutions is not physically realistic because of the formation of finite inversion/mixing layer in lower atmosphere that restrict the vertical pollutant diffusion. The surface based inversion in stable conditions influences most of the atmospheric phenomena such as the suppression of the dispersion of pollutants, formation of fog (Payra and Mohan, 2014), etc. Accordingly, the dispersion of pollutants in the presence of the surface based inversion is receiving significant attention in the research community (Ragothaman et al., 2001).
In an another approach, the concentration distribution c(x, y, z) of a pollutant released continuously from a point source in three-dimensional domain is given by the following relation (Huang, 1979;Buske et al., 2007;Irwin et al., 2007;Kumar and Sharan, 2009): where σ y is the diffusion parameter in crosswind (y) direction. Relation (1) considers the Gaussian plume in crosswind direction and a solution of the two-dimensional advectiondiffusion equation for crosswind integrated concentrations, C y (x, z), in vertical (z) and horizontal (x) directions. An advantage of this approach is that it required solving a twodimensional advection-diffusion equation for C y (x, z) that has comparably lesser complexity to find out a solution for the realistic profiles of wind speed and eddy diffusivity than a solution of the three-dimensional advection-diffusion equation. These types of dispersion models, derived from the analytical solution of two-dimensional advection-diffusion equation, have been used for air quality prediction (Kumar and Goyal, 2013); however, without any deposition of the pollutants on the ground surface. Recently, Tirabassi et al. (2008) and Kumar and Sharan (2014) derived the solutions of the two-dimensional advection-diffusion equation by considering the deposition on the ground surface. However, the mathematical techniques to solve the advection-diffusion equation in Tirabassi et al. (2008) and Kumar and Sharan (2014), and also in other numerical or analytical solutions are required to verify with an exact solution of this equation. In this study, an exact solution of the two-dimensional advection diffusion for crosswind integrated concentrations is derived for power law profiles of wind speed and vertical eddy diffusivity in the finite ABL. A dispersion model for a continuous point source release was proposed based on the derived analytical solution that consider the dry deposition on the ground surface. The effect of finite ABL on the dispersion of pollutant was accounted in the model through a reflecting boundary condition at top of the ABL. The model was evaluated with a depositing tracer observations from Hanford diffusion experiment in stable atmospheric conditions.

DEVELOPMENT OF THE DISPERSION MODEL
The crosswind integrated concentration C y (x, z) of a continuously releasing pollutant in the atmosphere in Eq. (1) is given by the two-dimensional steady state advectiondiffusion equation: where U(z) and K z (z) are the wind speed along the direction of horizontal mean wind (advection) and vertical eddy diffusivity profiles, respectively. In Eq.
(2), vertical component of the wind was neglected as its magnitude is much smaller than the horizontal wind components. The lateral component of wind was apparently included in the mean wind speed by rotating the horizontal x-axis in direction of the mean wind. In Eq.
(2), it was assumed that the scale of the horizontal wind speed is such as to render the longitudinal eddy diffusion negligible compared to the advection.
In order to derive an analytical solution of Eq.
(2), the profiles of the wind speed and eddy diffusivity are required to be specified. U(z) and K z (z) were approximated by the following empirical power law functions of vertical height above the ground surface where the power law exponents α and β in wind speed and vertical eddy diffusivity respectively are empirical parameters and U(z r ) and K z (z r ) are respectively the values of U(z) and K z (z) at a reference height z r . The values of power law exponents α and β in Eqs. (3) and (4) can empirically be estimated and depend on the atmospheric stability (Pasquill and Smith, 1983). The boundary conditions at the ground surface and top of the ABL specify the role of these boundaries and influence the characteristics of the pollutant dispersion in the ABL and need to be specified. A pollutant source can also be specified through a boundary condition. The physically realistic boundary conditions to solve Eq. (2) are imposed as follows: (i) A point source of strength Q was assumed to be located at a location (0, H s ), where H s is the effective height of that source. The source was continuously releasing the pollutant into the atmosphere. A point source can be represented in terms of the Dirac-delta function (δ) and by the following boundary condition at x = 0: (ii) The top of inversion/mixing layer located at height h acts like a potential blockade and limit the upward mixing of a pollutant. Thus, a perfectly reflecting flux boundary condition for the pollutant concentration was assumed at the mixing layer height and accordingly the vertical concentration flux vanishes at the height h is (iii) Incorporating the effects of removal or absorption of a contaminant at the ground surface can be defined through a deposition velocity. Thus, the vertical flux of a contaminant is proportional to the concentration at the ground surface and accordingly the boundary condition at z = 0 can be defined as (Calder, 1961): where V d is a proportionality constant with a unit of velocity and defined as the deposition velocity (V d ) of a pollutant at the ground surface. The sign of V d depends on whether the ground surface acts like a sink (positive) or a source (negative) (Koch, 1989) for a contaminant. If deposition velocity is zero (i.e., V d = 0) in Eq. (7), ground works as a perfectly reflecting surface to the pollutant. The deposition Eq. (7) is appropriate for small particles or gaseous pollutants because for heavy particles, one also needs to consider the gravitational settling velocity (w s ) in the advection diffusion Eq. (2). Then, the total deposition characterized by the quantity V d can be given by the sum of the turbulent transport along the concentration gradient and volume force effects causing uniform motion of the pollutant with a velocity w s in the vertical direction (Koch, 1989). However, consideration of the particle settling will increase the complexity of the problem to solve it analytically and thus not discussed in this study. Irrespective of the individual size and mass of the particles and other detailed mechanism, the deposition Eq. (7) considers their net effect at the ground surface.

Solution
An analytical solution of the two-dimensional steady state advection-diffusion equation (Eq. (2)) for crosswind integrated concentration C y (x, z) is derived using the separation-of-variables method. For power law profiles of the wind speed (Eq. (3)) and the vertical eddy diffusivity (Eq. (4)), Eq. (2) for C y (x, z) can be written as: Eq. (8) is a linear homogeneous parabolic partial differential equation (PDE) with linear homogeneous boundary conditions. Thus, a solution to the Eq. (8) was assumed in the following form of a product of two functions X(x) and Z(z): where X(x) and Z(z) are functions depending only on the independent variables x and z, respectively. Substituting of the sought solution (Eq. (9)) in Eq. (8), separate all the x's on one side of the equation and all the z's on the other side, and equating them to an arbitrary constant, it gives where λ 2 is called a separation constant. Eqs. (10) and (11) are linear ordinary differential equations (ODEs) of first and second order respectively and required to solve. The general solution of the X equation in Eq. (10) is: where A is an arbitrary constant. Substituting the Eq. (9) in boundary conditions (6) and (7), we find The second order ODE (Eq. (11)) with boundary conditions (13) and (14) form a regular Sturm-Liouville eigenvalue problem. The eigenvalues of this problem are real, positive, and discrete and there is a smallest eigenvalue. For each eigenvalue, an eigenfunction exists. The eigenfunctions corresponding to different eigenvalues are orthogonal with respect to the weighted inner-product with a weight function (a/b)z α (Zauderer, 1983). When V d ≠ 0, λ = 0 is not an eigenvalue. However, for V d = 0 in Eq. (14), zero will also be an eigenvalue and the corresponding eigenfunction is a constant, which may be taken as unity.
In order to estimate the non-zero eigenvalues and the corresponding eigenfunctions of problem (11), following transformation of the variables in Eq. (11) was applied These variables in Eqs. (15) and (16) transform the Eq. (11) to a following Bessel's equation: The general solution of Bessel's Eq. (17) yields Bessel functions of the first kind J ν and second kind J -ν and has the following form Substituting the solution (19) To estimate the arbitrary coefficients B 1 and B 2 in Eq. (20), the boundary condition at z = h (Eq. (13)) was applied. This boundary condition at top of the ABL with a recurrence relation J ν+1 (x) + J' ν (x) = νJ ν (x)/x (Abramowitz and Stegun, 1972) Similarly, the second boundary condition Eq. (14) at the ground surface z = 0 yields where  (.) is the gamma function. Eqs. (21) and (22) formed a set of two homogeneous linear equations and it will have infinitely many non-trivial solutions if its determinant vanishes and accordingly, one obtains and s n is given by Eq. (23) has non-trivial roots γ n ' and consequently the non-zero eigenvalues λ n can be determined from Eq. (18). For a γ n (eventually the eignevalue λ n ), substituting the expression of B 1 from Eq. (22) where the eigenfunction Z n (z) is given by Eq. (25) and the unknown coefficients A n (i.e., A 1 , A 2 , A 3 ,…… ) are required to be determine. The resulting solution derived in Eq. (26) characterizes the crosswind integrated concentration distribution C y (x, z) in horizontal and vertical directions.

Determination of A n in Eq. (26)
The unknown coefficients A n 's in solution (26) can be determines using the source condition (Eq. (5)) at x = 0. The concentration C y (0, z) in the solution (26) at x = 0 leads the source condition (5) to Since eigenfunctions Z n (z) are mutually orthogonal functions with respect to the weight function (a/b)z α (Abramowitz and Stegun, 1972), this property of orthogonally can be utilized to estimate A n . Multiply by an eigenfunction Z m (z), m ≥ 1 on both sides of Eq. (27) where pFq represents the generalized hypergeometric function.
The estimated coefficients A n from Eq. (28) where Z n (z) is given by Eq. (25) and Z n (H s ) is the value of the eigenfunction Z n (z) (Eq. (25)) at the source height at z = H s and R n is obtained by Eq. (29). This derived solution (Eq. (31)) is a closed form analytical expression for crosswind integrated concentrations released continuously from an elevated or ground level (by taking H s → 0 in the final solution Eq. (31)) point source in the atmospheric surface layer. The dispersion model incorporated the deposition of a pollutant at the ground surface. The concentration distribution in this expression (Eq. (31)) considered the effect of the finite ABL and varying wind speed and vertical eddy diffusivity.

Ground Level Releases
The effective source height H s can be assumed zero in Eq. (31) for a ground level point release. By taking H s → 0 in Eq. (31) and applying the properties of Bessel function for small arguments (Abramowitz and Stegun, 1972) where Z n (z) and R n are given by Eqs. (25) and (29), respectively. The expression (32) represents the concentration distribution released from a ground level point source.

Uniform Wind Speed and Eddy Diffusivity
For constant wind speed (i.e., U(z) = U ) and vertical eddy diffusivity (i.e., K z (z) = K ), the expression for crosswind integrated concentrations was derived by taking α = 0 (equivalently a = U ) and β = 0 (equivalently b = K ) in Eq. (31). For these α, β, a and b, Eq. (31) in which Eq. (23) for γ n 's transforms to This analytical solution (33) with Eq. (34) leads to a similar solution described by Fisher and Macqueen (1981) for the particulate matter transport from an elevated source in the atmosphere by assuming the zero gravitational settling velocity in their solution.

Without Deposition on the Ground Surface
When a pollutant is not depositing on the ground surface (i.e., V d = 0 m s -1 ), the zero vertical flux of the concentration in boundary condition (7) defines the ground surface to be act as a reflecting surface for the concentrations. By taking V d = 0 in Eq. (31), the dispersion model presented in section 2 reduces to that obtained by Lin and Hildemann (1996) for a perfectly reflecting ground surface for a contaminant. For a non-depositing tracer with zero deposition velocity at z = 0 in the lower boundary condition (7), zero will also be an eigenvalue and the corresponding eigenfunction is a constant and can be taken as unity. Thus, for V d = 0 m s -1 , the expression in Eq. (31) where γ n 's, deduced from Eq. (23), are the zeros of the following equation These expressions (Eqs. (35)-(36)) are same as obtained by Lin and Hildemann (1996) for crosswind integrated concentrations released from an elevated point source in the atmosphere. In this solution, both wind speed and eddy diffusivity are parameterized as the power law profiles of vertical height (Eq. (3) and Eq. (4)) and ground was considered to be a perfectly reflecting surface for a pollutant.

SENSITIVITY ANALYSIS OF THE MODEL WITH DEPOSITION VELOCITY
For an elevated release from an effective source height H s = 50 m, a sensitivity analysis of the present analytical dispersion model with different deposition velocities and with no-deposition (i.e., V d = 0 m s -1 ) on the ground surface was carried out. Meteorological parameters involved in power law profiles of the wind speed (Eq. (3)) and vertical eddy diffusivity (Eq. (4)) were arbitrarily taken as: α = 0.29, β = 0.45, z r = 10 m, U(z r ) = 2.1 m s -1 and K z (z r ) = 5 m 2 s -1 . The receptors were assumed to be at 2 m above the ground surface. Fig. 1 shows the direct influence of deposition velocity on the ground level normalized crosswind integrated concentrations released from a continuous elevated point source.
In Fig. 1, it is seen that up to a critical distance from source, normalized crosswind integrated ground level concentrations (glc) are increasing sharply, achieves maximum at a location, and beyond that it is decreasing with downwind distance. This behavior of the glc, released from an elevated point source of 50 m height, is obvious because regardless of the deposition velocity (V d ) at the ground surface, the maximum glc value and its occurrence location depend on the source height and the atmospheric stability. For a ground level point source, the glc is generally maximum at the source location and it decreases with the increasing downwind distances from the source. However, from the releases at higher vertical heights, the glc is smaller at near to the source and it increases with increasing downwind distance from the source and achieves maximum at a location and then further decreases (Fig. 1). Irrespective of the source height, the effect of the ground reflection (i.e., V d = 0) increases the ground level concentration; however, the absorption (i.e., V d > 0) at the ground surface reduces the glc. Thus, it is apparent that the partially absorption of pollutants at the ground surface reduces the ground level concentrations. The ground level crosswind integrated concentrations at different downwind distances from the source decrease with increasing deposition velocities (Fig. 1).

MODEL EVALUATION WITH A TRACER EXPERIMENT
In this study, Hanford diffusion experiment is utilized for evaluation of the present dispersion model for crosswind integrated concentrations. This tracer experiment was conducted at Hanford, USA during May-June, 1983. Two tracers (i) depositing tracer zinc sulfide (ZnS) and (ii) gaseous sulfur hexafluoride (SF 6 ) were released at 2 m height above the ground surface. Concentrations of both tracers were measured at five sampling arcs 100, 200, 800, 1600 and 3200 m downwind from the released point source (Doran et al., 1984;Doran and Horst, 1985). The receptors on each arc were placed at 1.5 m height above the ground surface. Doran and Horst (1985) reported the values of the meteorological and turbulence parameters like wind speed, friction velocity (u * ), inversion height (h) and the Obukhov length (L) and Table 1 presented these values for evaluation purpose. The values of deposition velocities V d at 800, 1600 and 3200 m arcs were also taken from Doran and Horst (1985). Doran and Horst (1985) computed these deposition velocities using a surface depletion method derived from a plume depletion model by Horst (1977).
The values of parameters α, β and K z (z r ) in power law profiles of wind speed and vertical eddy diffusivity in each trial are required for evaluation. These parameters for each trial are estimated from the observations of meteorological variables available in the Hanford diffusion experiment. To estimate the power law exponent α in the profile of wind speed in Eq. (3), a similar methodology used in Sharan and Kumar (2009) is utilized. The methodology is based on equating the values of first order derivatives of wind profile (Eq. (3)) and a profile given by the Monin-Obukhov similarity theory at a reference height z r . Consequently, it derived α in terms of the friction velocity (u * ), wind speed and stability function at a reference height z r .
The other power law exponent β in profile of vertical eddy diffusivity (Eq. (4)) was estimated by the Schmidt's conjugate law, β = 1 -α (Pasquill and Smith, 1983) and it is related to the estimated parameter α. The value of eddy diffusivity (K z (z r ) ) at a reference height z r in profile (3) was estimated using a flux-gradient relationship for power law wind profile (Eq. (3)) (Pasquill and Smith, 1983) in the surface layer in neutral to stable conditions. These parameters for each trial were computed from the meteorological observations in a similar manner as estimated in Sharan and Kumar (2009). Before evaluating the dispersion model  with real observations, a convergence analysis of the series in the proposed analytical solution (Eq. (31)) of advectiondiffusion equation was performed. The analysis was carried out by varying number of terms (N) in the series solution (31). For this analysis, it was observed that N = 100 is sufficient to reach the convergence of series in solution (31). The standard performance measures (Chang and Hanna, 2004) for evaluation of the air dispersion models were used in this study. These performance measures include the statistical indices: normalized mean square error (NMSE), factor of two (FAC2), fractional bias (FB), correlation coefficient (COR), and fractional variance (FS). COR and FAC2 would have unity and NMSE, FB, FS would have zero values for an idealized dispersion model (Chang and Hanna, 2004). Fig. 2(a) shows a scatter diagram between the ratios of predicted to observed crosswind integrated concentrations and the observed crosswind integrated concentrations by taking the deposition and without deposition on the ground surface in the proposed dispersion model. The crosswind integrated concentrations from model and observations were normalized by dividing them by the source strength Q. Fig. 2(a) shows that the simulated crosswind integrated concentrations from the proposed dispersion model with deposition on the ground surface are in good agreement with the observations. The proposed dispersion model with deposition simulated ~72% cases within a factor of two to the observed crosswind integrated concentrations. The NMSE, FB and FS between the computed and observed crosswind integrated concentrations are 0.20, -0.34 and -0.40 respectively. The modeled concentrations have a reasonably good correlation with the observations and it is also obvious from the value of correlation coefficient COR = 0.92. A negative fractional bias computed from the simulated and observed crosswind concentrations shows that proposed dispersion model is over-predicting from the observed concentrations. However, consideration a no-deposition condition in the model for a depositing tracer ZnS gives rise a larger over-prediction to the observations (FB = -1.05) in comparison to with deposition and all the points lie outside the factor of two lines ( Fig. 2(a)). Without considering the deposition in dispersion model, the NMSE and FS were found to 1.84 and -1.40, respectively. However, it was noted that the predicted concentrations with deposition and without deposition in the dispersion model have similar degree of correlation (COR = 0.92) with the observations. A quantilequantile (Q-Q) plot ( Fig. 2(b)) was also drawn by shorting the both observed and predicted concentrations in increasing order of their magnitudes (Venkatram, 1999). The Q-Q plot in Fig. (2b) also showed that the modeled crosswind integrated concentrations were close to the one-to-one line observations by taking deposition in comparison to without deposition on the ground surface.
The evaluation of the Hanford experiment in earlier studies, like in Doran and Horst (1985), is generally based on the four Gaussian plume-depletion models that required the vertical dispersion parameter, and also the parameterization of the constant form of the eddy diffusivity in form of the vertical dispersion parameter. However, the present approach is based on the variable forms of the vertical eddy diffusivity and wind speed in the atmosphere and does not required any specification of the dispersion parameters. However, the present solution required the computation of the hypergeometric function; though, with available computational resources in present days, it doesn't offer any difficulty in the computations.

CONCLUSIONS
A two-dimensional analytical dispersion model is proposed for the crosswind integrated concentrations released from a continuous point sources in the atmospheric surface layer with deposition to the ground surface. The proposed model considered the variability of the wind speed and eddy diffusivity as the power law functions of the vertical height above the ground surface. For these profiles of wind speed and eddy diffusivity, a closed form analytical solution of the resulting two-dimensional steady state advection-diffusion equation with deposition to the ground surface was derived using separation of variables technique. The model includes the effect of the finite ABL on dispersion. Several special cases of the proposed model were deduced for (i) a ground level point source release, (ii) uniform wind speed and eddy diffusivity, and (iii) without deposition on the ground surface. A sensitivity analysis of the model for ground level crosswind integrated concentrations with varying deposition velocities at the ground surface was performed. The derived solution of the advection-diffusion equation proposed in this study can be used to verify the mathematical techniques used in the previous or future studies.
For a depositing tracer ZnS, the proposed model was evaluated with the tracer observations obtained from Hanford diffusion experiment in stable conditions. The crosswind integrated concentrations computed from the present model with deposition are in good agreement with those observed. Model predicts ~72% cases within a factor of two. The statistical indices NMSE, FB and FS between the computed and observed concentrations are found to be respectively 0.20, -0.34 and -040 and model predicts a reasonably good correlation (COR = 0.92) with the observations. Present dispersion model was giving an over-predicting trend with the observed concentrations. However, considerations of a no-deposition condition in the dispersion model for a depositing tracer gives rise a larger over-prediction to the observations in comparison to with deposition. The derived exact solution of the advection-diffusion equation in this study can also be utilized to evaluate and error analysis of the numerical or other mathematical techniques for this equation.