## Abstract

For the dispersion of a pollutant released from a continuous source in the atmospheric boundary layer (ABL), a generalized analytical model describing the crosswind-integrated concentrations is presented. An analytical scheme is described to solve the resulting two-dimensional steady-state advection–diffusion equation for horizontal wind speed as a generalized function of vertical height above the ground and eddy diffusivity as a function of both downwind distance from the source and vertical height. Special cases of this model are deduced and an extensive analysis is carried out to compare the model with the known analytical models by taking the particular forms of wind speed and vertical eddy diffusivity.

The proposed model is evaluated with the observations obtained from Copenhagen diffusion experiments in unstable conditions and Hanford and Prairie Grass experiments in stable conditions. In evaluation of the model, a recently proposed formulation for the wind speed in the entire ABL is used. It is concluded that the present model is performing well with the observations and can be used to predict the short-range dispersion from a continuous release. Further, it is shown that the accurate parameterizations of wind speed and eddy diffusivity provide a significant improvement in the agreement between computed and observed concentrations.

## 1. Introduction

Dispersion of a continuous release of air pollutants in the atmospheric boundary layer (ABL) is an important problem owing to the complex nature of inhomogeneous sheared turbulence near a rough boundary, stability and meteorological conditions such as wind, temperature inversion and foggy atmosphere. When a plume is dispersed in the ABL, its shape, evolution and internal structure are determined by the interaction between the plume and the turbulent eddies that characterize the atmospheric motion. In spite of the recent advances in the theory and the direct numerical simulation of turbulent dispersion in the ABL, there remains a need of a fast flexible way to predict the transport of pollutants and other substances derived directly or indirectly from distributed surface sources in the atmosphere (Smith 2003).

The advection–diffusion equation (Seinfeld 1986) has long been used to describe the dispersion of air pollutants in the turbulent atmosphere. The use of analytical solutions of this equation was the first and remains the convenient way for modelling air pollution (Demuth 1978). The conventional air pollution model assumes wind speed and eddy diffusivity as constant throughout the entire ABL in deriving the solution of the advection–diffusion equation. In the solution obtained, the eddy diffusivity is parameterized in terms of dispersion parameters and downwind distance in analogy with the Gaussian distribution to describe the concentration close to the observations (Csanady 1973; Seinfeld 1986; Arya 1999). In stationary and horizontally homogeneous conditions, observational and theoretical studies show that the wind speed and eddy diffusivity vary with the vertical height above the ground (Stull 1988). Based on Taylor’s (1921) analysis and statistical theory, it is argued (Arya 1995; Demael & Carissimo 2008) that the eddy diffusivity depends on the downwind distance from the source (Mooney & Wilson 1993). Thus, over the years, the assumption of constant wind speed and eddy diffusivity has been relaxed to find the solution of the advection–diffusion equation.

Efforts have been made to obtain the analytical solutions of the two-dimensional steady-state advection–diffusion equation, describing the crosswind-integrated concentrations in the atmosphere, for particular forms of wind speed and vertical eddy diffusivity (Smith 1957; Van Ulden 1978; Huang 1979; Wilson 1982; Nokes *et al.* 1984; Turfus 1986; Lin & Hildemann 1996; Wortmann *et al.* 2005; Sharan & Modani 2006). The crosswind-integrated solution is often used to deduce the three-dimensional concentrations of a pollutant, assuming the Gaussian concentration distribution in crosswind direction, for environmental impact and assessment studies (Buske *et al.* 2007; Irwin *et al.* 2007). On the other hand, the crosswind-integrated concentration prohibits application to those problems where the concentration plume occupies a sufficiently deep region of the ABL to have been influenced by the height variation of the mean wind direction. In conditions of horizontally homogeneous turbulence in the ABL, most of the crosswind-integrated solutions (Smith 1957; Van Ulden 1978; Huang 1979; Lin & Hildemann 1996) are obtained by taking the wind speed and vertical eddy diffusivity as the power-law functions of vertical height above the ground. Turfus (1986) analysed the almost similar problem to that considered by Smith (1957) for dispersion from a line source in a reversed-flow layer by taking the power-law wind speed and later extended the analysis for a point source. Following the approach used by Smith (1982), Nokes *et al.* (1984) solved this equation for logarithmic wind speed and parabolic profile of eddy diffusivity and showed a significant improvement in the agreement between theory and observations. Wortmann *et al.* (2005) used the Laplace transform technique to solve the advection–diffusion equation for eddy diffusivity depending on both the downwind distance (*x*) from the source and vertical height (*z*) above the ground by dividing the *x*-domain into subdomains. In each subdomain, the averaged value of the eddy diffusivity and its derivative in the *x* variable is assumed. Recently, for near source dispersion, Sharan & Modani (2006) developed a model by taking the wind speed as a power-law profile and eddy diffusivity as a generalized function of downwind distance from the source.

In all these analytical models, the wind speed is either a power law or logarithmic profile of vertical height and similarly the eddy diffusivity has been assumed either a power law or a parabolic profile of *z* or a function of downwind distance from the source. But none of these provides a systematic approach to find the solution with the generalized functional forms of wind speed and eddy diffusivity. Even the solution given by Nokes *et al.* (1984) is limited to those forms of eddy diffusivity and wind speed that could accurately be expanded into the power series or polynomials and Wortmann *et al.* (2005) assumed constant averaged value of eddy diffusivity and its derivative in the *x* variable by performing a step-wise approximation. However, these profiles of wind speed and eddy diffusivity may not simulate the ABL processes realistically because of its complex turbulent structure in various stability conditions. More accurate and theoretically derived representations of the wind speed and eddy diffusivities in the ABL have been proposed in the literature (Monin & Obukhov 1954; Nieuwstadt 1984; Holtslag & Moeng 1991; Degrazia *et al*. 1997; Mooney & Wilson 1993; Gryning *et al.* 2007; and others) in various stability conditions. However, for these profiles of wind speed and eddy diffusivity in ABL, it is still not feasible to find a generalized analytical solution of the advection–diffusion equation that could effectively account for the influence of turbulent structure of the ABL on the dispersion of pollutants.

To overcome the limited efficacy of analytical dispersion models with the particular forms of wind speed and vertical eddy diffusivity to deal with the turbulent dispersion of a pollutant released from a continuous source in the ABL, a generalized analytical model describing the crosswind-integrated concentrations is presented. An analytical scheme is described to solve the resulting two-dimensional steady-state advection–diffusion equation for wind speed as a generalized function of vertical height above the ground and eddy diffusivity as a function of both downwind distance from the source and vertical height. This model preserves the beauty of analytical solution without compromising the accuracy of wind speed and eddy diffusivity to compute the crosswind-integrated concentrations.

## 2. Model description

### (a) Advection–diffusion equation

The two-dimensional steady-state advection–diffusion equation for the crosswind-integrated concentration *C*(*x*, *z*) of a non-reactive contaminant released from a source in the ABL based on *K*-theory is given by
2.1
where *U*(*z*) is the wind speed profile in *x*-direction and *K*_{x} and *K*_{z} are the longitudinal and vertical eddy diffusivities in the *x*- and *z*-directions, respectively. Equation (2.1) is derived by assuming unidirectional wind.

Under moderate to strong winds, the transport owing to longitudinal diffusion in equation (2.1) is neglected in comparison to advection and accordingly 2.2

### (b) Boundary conditions

Equation (2.2) is subject to the following boundary conditions:

(i) The ground is assumed to be a perfectly reflecting surface and accordingly diffusive flux vanishes close to the surface, i.e. 2.3

All the models in the literature assume the lower boundary at *z*_{b}→0. However, the consideration of *z*_{b}→0 as the lower boundary may be questionable as most of the parameterizations of wind speed and eddy diffusivity are not valid in the viscous sublayer. Most of the profiles of *K*_{z}(*x*,*z*) vanish as *z* tends to zero, implying that the boundary condition (2.3) becomes an identity satisfying for all vertical concentration gradients close to the surface.

It is evident that the most physical formulations of *U* and *K*_{z}, based on the Monin–Obukhov similarity theory (MOS), are applicable in the surface layer and no longer valid in the roughness sublayer. Micrometeorological measurements also indicate that the idealized surface layer of the similarity theory may not begin even at the top of the roughness elements of the average height *h*_{o}, which is an order of magnitude larger than *z*_{0} (where *z*_{0} represents the surface roughness length), but at 1.5–2 times *h*_{o} (Arya 1999), implying that the *z*_{b} should not be less than *z*_{0}. Accordingly, in the present study, we have taken *z*_{b}=*z*_{0}.

(ii) The pollutant is not able to penetrate through the top of the inversion/mixed layer located at height

*h*, i.e. 2.4(iii) The pollutant is released from an elevated point source of strength

*Q*located at the point (0,*H*_{s}) 2.5 where*H*_{s}is the stack height and*δ*is the Dirac-delta function.

### (c) Solution procedure

Equation (2.2) can be re-written as
2.6
Since the boundary conditions (equations (2.3) and (2.4)) are homogeneous, a basis of all linearly independent normalized vectors satisfying the eigenvalue problem
2.7a
2.7b
is given by
2.8
where *λ*_{n}=*n**π*/(*h*−*z*_{b}) (*n*=0,1,2,…) are the corresponding eigenvalues. The eigenvectors *Φ*_{n}’s in equation (2.8) are orthonormal. These eigenfunctions (equation (2.8)) in the form of the cosine waves arise in view of zero slop at *z*=*z*_{b} and *z*=*h*.

The analytical solution of equation (2.2) with the boundary conditions (equations (2.3) and (2.4)) has been obtained by taking the solution of the form
2.9
where *A*_{n}(*x*), *n*=0,1,2,… are the unknown coefficients of the series. To seek these coefficients *A*_{n}(*x*), we compel the solution (2.9) to satisfy the partial differential equation (2.6) and accordingly, we get
2.10
Multiplying both the sides by eigenfunctions {*Φ*_{m}(*z*)} (*m*=0,1,2,…) and integrating it with respect to *z* from *z*_{b} to *h*, we obtain a following set of first-order ordinary differential equations (ODEs)
2.11
This can be re-written in the matrix notation
2.12
in which *A*(*x*)=[*A*_{0}(*x*),*A*_{1}(*x*),*A*_{2}(*x*),*A*_{3}(*x*),…]^{t} (here the superscript t denotes the transpose) is the column vector of unknown coefficients and
2.13
where the elements of matrices *B* and *E* are, respectively,
2.14a
2.14b
In equations (2.13) and (2.14), matrices *B* and *E* appear to be rectangular matrices for different upper integer bounds of *m* and *n*, which seems to require the generalized inverse of a rectangular matrix *B* (Greville 1959; Harwood *et al.* 1970). Without the loss of generality, to avoid the complicated procedure involved in finding the generalized inverse of a rectangular matrix *B*, we vary the *m* and *n* up to the same upper integer value *N* for which both *B* and *E* become the square matrices of order *N*. *B* is a symmetric and positive definite matrix with non-zero positive diagonal elements (electronic supplementary material 1) implying that the matrix *B* is invertible.

In equation (2.12), *F* is an *N*×*N* matrix of the real-valued continuous functions of real variable *x*. The system of first-order ordinary linear differential equations (2.12) will have a unique solution (Martin 1967; Zhu 1992)
2.15
provided the matrix *F*(*x*) satisfies the condition
2.16
where the column matrix *A*(0) is the value of the matrix *A*(*x*) at *x*=0. The functional forms of *K*_{z}(*x*,*z*) used in the dispersion modelling for which the corresponding matrix *F* satisfies the condition (2.16) are deduced in appendix A.

To determine *A*(0), solution (2.9) is substituted into the source condition (equation (2.5))
2.17
Multiplying both sides of equation (2.17) by eigenfunction *Φ*_{m}(*z*) and integrating it with respect to *z* from *z*_{b} to *h*, one obtains
2.18
or in matrix notation
2.19
where *G*=*Q**Φ*_{m}(*H*_{s}) is a *N*×1 column matrix and *B* is given by equation (2.14a).

### (d) Computation of exponential of matrix *F*(*x*)

For all the functional forms of *K*_{z}(*x*,*z*) derived in the appendix A, the matrix *F*(*x*) in equation (2.15) can be expressed as
2.20
where *M* is a square matrix with elements independent of the variable *x*, henceforth refers to a constant matrix. The functional form *f*(*x*) arises corresponding to the forms of *K*_{z} as *K*_{z}(*x*,*z*)=*f*(*x*) and *K*_{z}(*x*,*z*)=*f*(*x*)*g*(*z*). The matrix *M* can be transformed into the following Jordan canonical form (Horn & Johnson 1985):
2.21
where *J* is the Jordan matrix of all eigenvalues of *M* and *P* is the corresponding matrix of all linearly independent eigenvectors. If all the eigenvalues of *M* are distinct, *J* will be a diagonal matrix (Horn & Johnson 1985). Thus, the integration of equation (2.20) can be written as
2.22
where
2.23
is a diagonal matrix. By taking the exponential of both sides of equation (2.22), one gets
2.24
in which the exponential of diagonal matrix *D* can be obtained by simply taking the exponential of the diagonal elements of *D*.

Once the diagonal elements of the matrix *D* are known, the expression on the left-hand side of equation (2.24) becomes known and in turn the coefficients {*A*_{n}(*x*)} will be evaluated from equation (2.15). Consequently, these coefficients’ *A*_{n}’s and the eigenfunctions in the series (equation (2.9)) allow to compute the concentration *C*(*x*,*z*).

For a generalized *F*(*x*) in equation (2.15), satisfying the condition (2.16), it is not obvious to find the variable eigenvalues and the corresponding eigenvector matrix for the matrix . However, from the computational point of view, we can calculate these eigenvalues and the corresponding eigenvector matrix for any preset value of *x* for which the coefficients of the matrix become constants. Consequently, the exponential of the matrix in equation (2.15) can be calculated for any preset value of *x*.

We have derived the solution of equation (2.2) by assuming the lower boundary at *z*=*z*_{b}, avoiding the problem associated with the vanishing of *K*_{z} at the ground *z*=0. For the dispersion in the ABL, *z*_{b} may be taken as *z*_{0}, a roughness length scale. The solution with the lower boundary condition at the ground *z*=0 can be readily reduced by simply taking the limit as *z*_{b}→0. The solution for a ground-level release is obtained when the effective stack height *H*_{s} tends to zero in equation (2.5). The form of the solution corresponding to a ground-level release is similar to that of elevated release except the term *Φ*_{m}(*H*_{s}) tends to one in equation (2.19) and consequently the column matrix *G* in equation (2.19) becomes *G* = *Q*. We hasten to point out that the analogy of the approach used here is similar to the Galerkin technique in which the basis functions are chosen satisfying the boundary conditions.

## 3. Comparison with existing models

In order to check the accuracy of the employed analytical solution technique given in §2, the solution of equation (2.2) given by equation (2.9) is reduced into some particular cases. However, it is not always feasible to reduce in the form of an explicit expression from the solution given in §2*c* for all the parameterizations of *U*(*z*) and *K*_{z}(*x*,*z*). Thus, the proposed model is compared with other analytical models obtained by taking the particular forms of *U*(*z*) and *K*_{z}(*x*,*z*). As most of the analytical models in the literature are derived by considering the lower boundary at the ground, we have taken *z*_{b}=0 for the comparison purposes.

Most of these analytical models can be broadly classified first into two groups based on the wind field as (i) uniform wind speed and (ii) variable wind speed—a function of *z*. In each of these groups, the model differs with the choice of the eddy diffusivity *K*_{z}(*x*,*z*).

### (a) Uniform wind speed

#### (i) Constant eddy diffusivity *K*_{z}(*x*,*z*)=*K*

For constant *U* and *K*, the matrix *B* becomes a scalar matrix with the scalar *U* and *E* is a diagonal matrix with the elements
3.1
Thus, the matrix *F* in equation (2.15) will be a diagonal matrix with diagonal elements
3.2
In equation (2.15), the diagonalization of matrix *F* instantly calculates the integration of *F* by just integrating the diagonal elements *f*_{nn} and then taking their exponential leading to
3.3
The coefficients *A*_{n}(0) in equation (2.19) are given by
3.4
Thus, by using *λ*_{n} and *Φ*_{n} in equation (2.8) with *z*_{b}→0, the solution of equation (2.2) for constant *U* and *K*_{z} transforms equation (2.9) to
3.5
This is the classical solution (Seinfeld 1986) which leads to the Gaussian plume model for an elevated source with reflections from two parallel boundaries at *z*=0 and *z*=*h*.

#### (ii) *K*_{z}(*x*,*z*) as a function of downwind distance

For the treatment of the near source dispersion, eddy diffusivity is parameterized in terms of downwind distance *x* (Arya 1995; Sharan *et al.* 1996). In a similar manner to case §3*a*(i), for a uniform wind speed and *K*_{z}(*x*,*z*)=*f*(*x*), a function of downwind distance, the solution described in §2*c* is reduced to the following form:
3.6
This solution (3.6) is the same as that obtained by Sharan & Gupta (2002) for uniform wind speed and eddy diffusivity as a generalized function of *x*. If the eddy diffusivity in equation (3.6) is taken as constant, equation (3.6) reduces to equation (3.5). This model has been evaluated (Sharan & Gupta 2002) with the Copenhagen tracer observations by taking *f*(*x*) as a linear function of downwind distance.

### (b) Variable wind speed as a function of vertical height above the ground

#### (i) Power-law wind speed and eddy diffusivity in terms of downwind distance

For predicting the crosswind-integrated concentration, wind speed and eddy diffusivity are parameterized as (Pasquill & Smith 1983)
3.7
3.8
where *U*(*z*_{r}) is the wind speed at a reference height *z*_{r} and *α* is the power-law exponent that depends on the atmospheric stability and can be parameterized by MOS theory (Huang 1979). The analytical solution of equation (2.2) with boundary conditions (2.3)–(2.5) (assuming *z*_{b}=0) with the parameterization of wind speed (equation (3.7)) and diffusivity (equation (3.8)) is obtained by the method of eigenfunction expansion (Sharan & Modani 2006) and is given by
3.9
where *μ*=1/(*α*+2), *J*_{−μ} is the Bessel function of the first kind and *β*_{n}’s are the roots of
3.10
For computational purpose, (Sharan & Modani 2006), in which (Arya 1995), where *σ*_{w} refers to the standard deviation of the vertical velocity component and is the average wind speed.

For the comparison, the effective stack height of an elevated source is taken 50 m and the mixing height is assumed to be 500 m. We have taken the parameters *α*=0.29, *U*(*z*_{r})=2.1 ms^{−1}, *z*_{r}=10 m in the power-law profile of wind speed and *σ*_{w}=0.83 (ms^{−1}), (ms^{−1}) ( is computed by taking the average of power-law wind profile (equation (3.7)) in ABL) in *f*(*x*).

Figure 1*a* compares the concentrations computed by the present technique to that from the closed form analytical solution (3.9) for a ground-level receptor at *z*=0.02 m and at plume centre line at *z*=50 m with increasing *x*. This shows that the concentrations computed from both the techniques overlap each other for all *x*. For *z*=0.02 m, the maximum absolute error is found to be 0.5 per cent occurring near to the source. However, this error is attributed to computational errors in both of the solutions (2.9) and (3.9).

#### (ii) Both wind speed and eddy diffusivity in terms of power-law functions of vertical height

Lin & Hildemann (1996) (after this LH96) proposed a scheme to solve equation (2.2) by using Green’s function theory for both wind speed and eddy diffusivity in terms of the power-law functions of vertical height above the ground. In their solution, the power-law profile for *U*(*z*) is given by equation (3.7) and eddy diffusivity is taken in the form
3.11
where *K*_{z}(*z*_{r}) is the eddy diffusivity at the reference height *z*_{r}. The analytical solution of equation (2.2) for these forms of *U*(*z*) and *K*_{z}(*x*,*z*) is (LH96)
3.12
with
3.13
where *γ*_{n}’s are the zeros of the following equation:
3.14
In the power-law profile of *K*_{z}(*x*,*z*) in equation (3.11), we have taken the parameters *β*=0.45, *K*_{z}(*z*_{r})=5.0 m^{2} s^{−1} and *α*, *U*(*z*_{r}) and *z*_{r} are given in the previous case (§3*b*(i)), figure 1*b* compares the computed concentrations from equation (2.9) to that with LH96 (equation (3.12)) and show that both almost overlap each other for all *x* except very near to the source at *z*=0.02 m, where the maximum absolute error is found to be less than 0.2 per cent.

#### (iii) Wind speed in terms of power-law functions of vertical height and K_{z}(*x*,*z*) is an explicit function of both downwind distance and vertical height

For near source dispersion, eddy diffusivity is also a function of downwind distance *x*. Mooney & Wilson (1993) have proposed a modified functional form of eddy diffusivity in terms of both the downwind distance *x* and the height *z* above the ground to deal with the near source dispersion. The modified form of *K*_{z}(*x*,*z*) is given in the form (Mooney & Wilson 1993)
3.15
in which *g*(*z*) is any form of eddy diffusivity depending on the *z* and *f*(*x*) is the correction to *g*(*z*) for near source dispersion and is a dimensionless integrable function of *x*. A power-law parameterization of *g*(*z*), given by equation (3.11), is taken in equation (3.15). Using the power-law profile of *U*(*z*) (equation (3.7)) and the modified parameterization of *K*_{z}(*x*, *z*) (equation (3.15)), a closed form analytical solution of equation (2.2) with boundary conditions (2.3)–(2.5) (assuming *z*_{b}=0) has been obtained by the method of separation-of-variables as (Sharan & Kumar 2009)
3.16
where *μ* is given by equation (3.13) and *γ*_{n}’s are the zeros of equation (3.14).

This solution is the generalization of SM06 and LH96 for near source dispersion. By taking *β*=0 and *b*=1.0 in equation (3.11), the solution (3.16) is reduced to SM06 (equation (3.9)). Similarly, the solution LH96 can be obtained by taking the function *f*(*x*)=1.0 (as a constant function).

For comparing the solutions (equation (2.9)) and equation (3.16), *f*(*x*) is taken a dimensionless function of *x* from Mooney & Wilson (1993). Figure 1*c* reveals that the computed concentrations from the proposed technique almost exactly match those obtained analytically from equation (3.16) for all *x* except near to the source. The maximum error is found to be less than 5 per cent (figure 1*c*) for near to the source at *z*=0.02 m. This error is attributed to the computational errors including round-off involved owing to the number of terms taken in the series solution calculating the concentration in each model.

#### (iv) Wind speed and eddy diffusivity are generalized functions of vertical height above the ground

In the above cases, wind speed has been parameterized as a power-law profile of vertical height above the ground. However, in the atmospheric surface layer, a more realistic profile of wind speed can be obtained by using the MOS or the local similarity theory (Nieuwstadt 1984). Also, *K*_{z}(*x*,*z*) can have a cubic polynomial of *z* (O’Brien 1970) in unstable ABL and a quadratic polynomial in stable ABL (Nieuwstadt 1984), vanishing at the ground as well as at the top of the mixing layer (Nokes *et al.* 1984).

In order to compare the solution (2.9) of equation (2.2) for a generalized profile of wind speed and eddy diffusivity, we have used the solution procedure proposed by Nokes *et al.* (1984) to the problem of turbulent flow in a long wide uniform channel. The procedure of separation-of-variables in equation (2.2) leads to a system of two ODEs in each of the variables *x* and *z*. The resulting second-order ODE in *z* is solved as a Sturm–Liouville problem using the power-series technique to determine the eigenvalues and associated eigenfunctions for the logarithmic profile of wind speed and parabolic eddy diffusivity. A slightly changed profile of *U*(*z*) from Nokes *et al.* (1984) to avoid the negative values of wind speed near to the surface is taken. Accordingly, we have assumed the following forms of logarithmic wind speed and parabolic eddy diffusivity:
3.17
3.18
where is the depth-averaged velocity, is the dimensionless parameter, is the diffusivity coefficient and *z*_{0m} is chosen such that *U*(*z*) becomes zero at *z*=0.

For comparison, we have taken the parameter’s values in profiles (3.17) and (3.18) from Nokes *et al.* (1984). Figure 1*d* demonstrates that the results computed using the technique described in this paper are in good agreement to those obtained from the one based on the power-series method adopted by Nokes *et al.* (1984) for *H*_{s}=50 m and *h*=500 m at *z*=0.02 and *z*=50 m. The maximum absolute error for ground-level receptor is found to be 2 per cent.

## 4. Results and discussion

The model given here is general in the sense that it can use any form of parameterization of wind speed *U*(*z*) and eddy diffusivity *K*_{z}(*x*,*z*) as an explicit function of variables *x* and *z*. The comparison results demonstrate that the crosswind-integrated concentrations computed by the proposed technique are in very good agreement with those obtained from the analytical solutions corresponding to the particular functional forms of *U*(*z*) and *K*_{z}(*x*,*z*). The computed concentrations from the proposed technique almost exactly match for all *x* except near to the source.

An analysis for the stability of the solution (2.9) is given in the electronic supplementary material 2. The mathematical analysis shows that the present solution is asymptotically stable. Also, based on the computations (electronic supplementary material 3), it is observed that *N*=250 terms are adequate to obtain a converged solution (equation (2.9)).

### (a) Model evaluation

The performance of the present model has been analysed from the Copenhagen diffusion experiment in unstable conditions and the Hanford and Prairie Grass experiments in stable conditions. We have used the standard statistical performance measures (i) normalized mean square error (NMSE), (ii) fractional bias (FB), (iii) correlation coefficient (COR), (iv) fractional variance (FS), and (v) factor of two (FAC2) that characterize the agreement between model prediction and observations (Sharan & Kumar 2009). A perfect model would have the idealized values as NMSE=FB=FS=0 and COR=FAC2=1.

#### (i) Copenhagen diffusion experiment in convective conditions

A tracer experiment in unstable and neutral conditions was conducted in the north Copenhagen area (Gryning & Lyck 1984; Gryning *et al.* 1987). In this experiment, a tracer SF_{6} was released without buoyancy at a tower height of 115 m and samples were collected on different arcs at 2 m above the ground surface. A detail description of meteorological and tracer data is provided by Gryning & Lyck (1984) and Gryning *et al.* (1987). The values of parameters such as *h*, friction velocity (*u*_{*}), Obukhov length (*L*), *σ*_{v}, *σ*_{w} and *U* observed during the experiment are taken from Gryning *et al.* (1987). The site was primarily residential with a roughness length of 0.6 m and thus the terrain was considered as an urban terrain in the computation.

*Parameterization of U(z) in unstable condition*

Normally, *U*(*z*) is parameterized in the dispersion model using an empirical formulation in terms of power-law function of *z*. The wind speed profile *U*(*z*) is also parameterized based on surface layer theory and Monin–Obukhov scaling. However, the profile based on the surface layer theory is valid only within the surface layer. This profile has been extended in the entire ABL by assuming a constant profile of wind speed above the surface layer. Recently, Gryning *et al.* (2007) have proposed a formulation for the wind speed profile in the entire ABL by extending the surface layer profile of wind speed over homogeneous terrain beyond the surface layer. In this profile, the local length scale in wind speed is composed of three component length scales. In the surface layer, the first length scale (*L*_{SL}) is taken to increase linearly with height with a stability correction following MOS theory. Above the surface layer, the second length scale (*L*_{MBL}) becomes independent of height but not stability, and at the top of the boundary layer, the third length scale is assumed to be negligible.

For unstable atmospheric conditions, the wind speed profile is given as (Gryning *et al.* 2007)
4.1
where
4.2
in which *ξ*=(1−12 *z*/*L*)^{1/3} and *L*_{MBL} is given by
4.3
where *f*_{c} is the Coriolis parameter, *k*=0.4 is von Karman constant and *A* and *B* are the resistance law functions depending on the dimensionless stability parameter *μ*=*u*_{*}/*f*_{c}*L* in the steady-state ABL, or *h*/*L* in the evolving ABL. For the neutral atmosphere, the values of *A*≈4.9 and *B*≈1.9 are rather well established. However, in stable and unstable conditions, the stability dependence of these functions *A* and *B* is a matter of discussion (Zilitinkevich & Esau 2005). Owing to the ambiguity in the formulation of the functions *A* and *B*, Gryning *et al.* (2007) have also proposed an empirical formula for the length scale *L*_{MBL} in both unstable and stable conditions as
4.4
In the present study, *L*_{MBL} is taken from equation (4.4) to compute the wind profile from equation (4.1).

*Parameterization of K _{z}(x,z) in unstable conditions*

In the dispersion models, *K*_{z} is parameterized as a function of the height above the ground, *z*, only. For the near source dispersion, *K*_{z} has also been parameterized as a function of downwind distance *x* (Arya 1995; Sharan & Modani 2006). Mooney & Wilson (1993) have used a modified form of *K*_{z} as a function of both variables *x* and *z* as
4.5
where *K*′_{z} is the functional form of eddy diffusivity depending on *z* only and *L*_{1}=*U*(*H*_{s})*τ*(*H*_{s}) is the along-wind length scale across which the diffusivity achieves its ‘far-field’ value and *τ* is the Lagrangian time scale at the source height.

The following form of *K*′_{z}(*z*) for unstable conditions as a function of *z* derived by Degrazia *et al.* (1997) in terms of the convective scaling parameters is used in the present study in formulation (4.5)
4.6
where *w*_{*} is the convective velocity scale and computed from *u*_{*} using the relation *w*_{*}=*u*_{*}(−*h*/(*k**L*))^{1/3} (Seinfeld 1986). The Lagrangian time scale *τ* at the source height *z*=*H*_{s} is parameterized as (Mooney & Wilson 1993)
4.7
where *σ*_{w} is the vertical turbulent intensity.

Notice that in the model, effectively, *K*_{z}(*x*,*z*) has been parameterized as a function of both *x* and *z*. However, in order to analyse the contribution arising from the dependency of *K*_{z}(*x*,*z*) on *x* in equation (4.5), we have also parameterized *K*_{z}(*x*,*z*)=*K*′_{z}(*z*) as a function of *z* only by equation (4.6). Figure 2*a* shows that the normalized crosswind-integrated concentrations computed from the present model in unstable conditions are in good agreement with the observations. All the points are predicted within a factor of two. In addition, it reveals that the predicted concentrations for both the parameterizations of *K*_{z}—depending on either both *x* and *z* or *z* only—are comparable to each other. In fact, the predicted concentrations from *K*_{z}(*x*,*z*) appear to be relatively closer to the observations than those obtained with *K*_{z} as a function of *z* only. It is noticed that the highest concentrations correspond to the shortest downstream distances, where the *x*-dependency of the *K*_{z} is expected to be significant. The extent of the dependency of the parameterization of *K*_{z} on *x* can be analysed by the contribution of the term *e*^{−x/L1} (equation (4.5)). When *H*_{s}→0 corresponding to a ground-level release, it leads to *L*_{1}→0 and the contribution of the term *e*^{−x/L1} becomes negligible. It implies that *K*_{z}(*x*,*z*) reduces to *K*′_{z}(*z*) and *K*_{z}(*x*,*z*) will be a function of *z* only. For *H*_{s}≠0, the value of *L*_{1} increases as *H*_{s} increases and thus the contribution of the term *e*^{−x/L1} will appear in the parameterization of *K*_{z}(*x*,*z*) in terms of both *x* and *z* for an elevated source. For *H*_{s}=115 m in this experiment, the difference between the results for the parameterizations of *K*_{z}(*x*,*z*) and *K*′_{z}(*z*) is clearly seen (figure 2 and table 1) and it shows that the dependency of *x* in the *K*_{z} gives a relatively better result. The unpaired analysis, *Q*–*Q* plot (figure 2*b*) drawn by arranging both the observed and the predicted concentrations in increasing order of their magnitudes (Venkataram 1999), shows that the predicted concentrations are close to the one-to-one line observations for both parameterizations of *K*_{z}(*x*, *z*). This fact is re-affirmed from the statistical measures given in table 1. The overall statistical measures indicate that the present model has a good correlation with the observations.

We have also compared the results from the present model with a model obtained by the generalized integral transform technique (GILTT) given by Wortmann *et al.* (2005). For this purpose, we have taken the same parameterizations of the wind speed and eddy diffusivity as used by Wortmann *et al.* (2005). Statistical indices in table 1 indicate that the results from the present model are comparable with those obtained using GILTT.

#### (ii) Hanford tracer experiment in stable conditions

In stable conditions, a diffusion experiment was conducted in May–June, 1983, at the Hanford site in eastern Washington (46^{°}34 N, 119^{°}36 W) USA (Doran & Horst 1985). The terrain at the Hanford Washington site is flat with a displacement height of approximately 1.4 m (Doran & Horst 1985). In the SF_{6} tracer experiment, six test runs were conducted in 1983 and the concentrations of SF_{6} were measured at five sampling arcs 100, 200, 800, 1600 and 3200 m downwind from a fixed release point. The tracer was released at 2 m and was sampled at 1.5 m height above the ground. Values of the meteorological parameters *h*, *u*_{*} and *L* in the evaluation of the model are taken from Doran *et al.* (1984). Values of surface roughness length *z*_{0} for each run are adopted from the re-analysis of the data by Kim & Larson (2001).

*Parameterization of U(z) in stable conditions*

For stable conditions, the extended profile of wind speed into entire ABL is given as (Gryning *et al.* 2007)
4.8
where *b*=4.7 and *L*_{MBL} in stable conditions is
4.9
As mentioned earlier, owing to the uncertainty involved in the determination of resistance law functions *A* and *B* as a function of *μ* in stable conditions, we have adopted the empirical relation (4.4) for *L*_{MBL} to calculate the wind profile from equation (4.8) in the present study.

*Parameterization of K_{z}(x,z) in stable conditions*

For dispersion in stable conditions, *K*_{z}(*x*,*z*) is parameterized by the same expression given in equation (4.5). However, the following form of *K*′_{z}(*z*), derived by Degrazia & Moraes (1992) from the local similarity theory and the statistical diffusion theory, is taken:
4.10
where *Λ* is the local Obukhov length, defined by *Λ*=*L*(1−*z*/*h*)^{(3α1/2)−α2}, in which *α*_{1} and *α*_{2} are constant parameters. We have used *α*_{1}=3/2 and *α*_{2}=1, derived by Nieuwstadt (1984) from the second-order closure model for stationary turbulence with constant cooling rate. In stable conditions, *σ*_{w} in equation (4.7) is computed by the following expression given by Nieuwstadt (1984):
4.11
In order to analyse the contribution arising from the dependency of *K*_{z}(*x*,*z*) on *x* in equation (4.5), we have also parameterized *K*_{z}(*x*,*z*)=*K*′_{z}(*z*) as a function of *z* only by equation (4.10). Figure 3*a* shows that the normalized crosswind-integrated concentrations computed from the present model in stable conditions are in good agreement with the observations and it predicts 73 per cent cases in a factor of two to observations. In addition, it reveals that the predicted concentrations for both the parameterizations of *K*_{z}-depending on either both *x* and *z* or *z* only—are comparable to each other. This is justifiable as both the parameterizations become identical because of the negligible contribution of the term *e*^{−x/L1} in equation (4.5) for ground-level releases. The unpaired analysis, *Q*–*Q* plot in figure 3*b* shows that the predicted concentrations are close to the one-to-one line observations for both of the parameterizations of *K*_{z}(*x*,*z*). This is re-affirmed from the values of statistical measures given in table 1.

We have also compared the results from the present model with a model for near source dispersion proposed by Sharan & Modani (2006) (SM06). The overall statistical measures in table 1 indicate that the present model predicts reasonably better than SM06. The NMSE of the present model is 0.20 whereas for SM06 it is 0.52. The COR shows that the present model has a good correlation with the observations. It is worth noting that the model leads to negative FB, i.e. a general overestimation of the concentration data is pointed out.

#### (iii) Prairie Grass experiment

The Prairie Grass experiment (Barad 1958) was conducted in O’Neill, Nebraska, during the summer of 1956 and is a standard and widely used dataset of dispersion observations used for the evaluation of dispersion models for ground-level releases of continuous plume over flat terrain. In all 68 test runs of this experiment, a tracer sulphur dioxide (SO_{2}) was released without buoyancy at a height of 0.46 m (except four runs# 65–68 in which the release height was 1.5 m) and the samples of the concentrations of SO_{2} were measured at 1.5 m above the ground surface on five sampling arcs 50, 100, 200, 400 and 800 m downwind from the source. We have used the digital version of the dataset (http://www.dmu.dk/International/Air/Models/Background/ExcelPrairie.htm) for normalized crosswind-integrated concentrations and meteorological variables. In this dataset, the meteorological parameters *L*, *u*_{*}, *w*_{*} and *h* are given and the surface roughness length *z*_{0} was reported as 0.006 m.

In the evaluation, the computations are carried out with this dataset containing the runs in all stability conditions from very unstable to neutral and very stable. However, we report the results for the runs in stable conditions only. In this stability, we have also abandoned the runs of very stable conditions (*h*/*L*>10) in which the turbulence becomes very weak, sporadic, intermittent and no longer continuous in time and space (Holtslag & Nieuwstadt 1986). In computations, the release and sampling heights are taken as 0.46 and 1.5 m, respectively. We have used the extended wind profile (equation (4.8)) into the entire ABL and the vertical eddy diffusivity profile is taken as the same as that used in the evaluation of the Hanford diffusion experiment in stable conditions. The values of *σ*_{w} in equation (4.7) are computed from the expression (4.11).

The scatter diagram between the ratios of predicted to observed crosswind-integrated concentrations and the observed concentrations is shown in figure 4*a*. In order to analyse the contribution arising from the dependency of *K*_{z}(*x*,*z*) on *x* in equation (4.5), we have also parameterized *K*_{z}(*x*,*z*)=*K*′_{z}(*z*) as a function of *z* only (equation (4.10)). Figure 4*a* shows that the predicted concentrations from the present model are in good agreement with the observations. The predicted concentrations for both the parameterizations of *K*_{z}—depending on either both *x* and *z* or *z* only—are comparable to each other. This is expected due to the negligible contribution of the term *e*^{−x/L1} in equation (4.5) for a ground-level release. For both of the parameterizations of *K*_{z}(*x*,*z*), the model predicts 95 per cent cases in the factor of two to observations and the NMSE and FB between the computed and the observed concentrations are, respectively, 0.06 and −0.16 (table 1). The *Q*–*Q* plot (figure 4*b*) reconfirms that the concentrations predicted from the present model are close to the one-to-one line for both parameterizations of *K*_{z}(*x*,*z*).

The results discussed here in both unstable and stable conditions are for *U*(*z*) based on the empirical parameterization of *L*_{MBL} in equation (4.4) which was obtained from different datasets and is observed better for stable than unstable conditions (Gryning *et al.* 2007). Owing to the inadequate meteorological data in the Copenhagen, Hanford and Prairie Grass experiments, we could not find the appropriate resistance law functions *A* and *B* to compute *L*_{MBL} required in the computation of *U*(*z*). However, in unstable conditions for the Copenhagen diffusion experiment, we tried to use the empirical relations of functions *A* and *B* proposed by Arya (1975) based on the data from the two sites (O’Neill in Nebraska, USA, and Hay in New South Wales, Australia). For *K*_{z}(*x*,*z*)=*K*′_{z}(*z*) in unstable conditions, the improved parameterization of *L*_{MBL} in the wind profile based on these resistance law functions reveals that the NMSE and FB reduce from 0.086 to 0.05 and 0.12 to 0.036, respectively. It also gives a better correlation (COR=0.92) with the observations in comparison to the COR=0.90 for empirical relation of *L*_{MBL} (equation (4.4)). It shows that a realistic parameterization of wind profile can lead to a significant improvement in the predication of concentrations of pollutants in the ABL.

Here, in comparison with measurements, it is noted that the model tends to underestimate (FB>0) for unstable conditions, but to overestimate (FB<0) for stable conditions. The underpredicting (overpredicting) trends of the concentrations from the model in unstable (stable) conditions may be attributed to a slight overestimation (underestimation) of the computed mean wind speeds from the formulation of Gryning *et al.* (2007) in comparison to those observed.

Here, the model is evaluated in stable conditions with the observations corresponding to the releases from a source near to the ground. On the other hand, the model is proposed here for both elevated and ground-level sources. Thus, the model needs to be evaluated further once the diffusion data become available for an elevated source in stable conditions.

### (b) Three-dimensional concentration at a specific location in the atmosphere

The concentrations obtained from solution (2.9) are two-dimensional and the crosswind dispersion term in advection–diffusion equation (2.2) is not explicitly treated. However, if one is interested in the impact or exposure studies of the concentrations at a given location (*x*, *y*, *z*) in the atmosphere rather than the concentration integrated along the crosswind direction *y*, the lateral diffusion term needs to be included to calculate the three-dimensional concentration *c*(*x*,*y*,*z*). Very limited analytical solutions are available by taking the lateral diffusion term in the advection–diffusion equation even for the specific forms of wind speed and eddy diffusivities (Yeh & Huang 1975; Huang 1979). However, based on the analysis of Prairie Grass and some other historical tracer experiments of atmospheric dispersion, Irwin *et al.* (2007) concluded that the ground-level crosswind concentration profile of a dispersing plume on average is well characterized as having a Gaussian shape, which is also predicted by all atmospheric transport and diffusion models, regardless of their sophistication. Thus, by assuming the Gaussian concentration distribution in crosswind direction (Huang 1979; Buske *et al.* 2007; Irwin *et al.* 2007), the steady-state concentration of a pollutant released from a point source at (0, 0, *H*_{s}) in a three-dimensional domain can be described as
4.12
where *C*(*x*,*z*) is the crosswind-integrated concentration and can be obtained from equation (2.9) for any generalized profiles of *U*(*z*) and *K*_{z}(*x*,*z*) and *σ*_{y} is the standard deviation of the concentration distribution in crosswind direction (Huang 1979; Pasquill & Smith 1983).

### (c) Issues and limitations

So far, an analytical dispersion model describing the crosswind-integrated concentrations is presented and evaluated with the observations and its application to the three-dimensional dimension is also described. The present work is subject to the underlying assumptions and limitations (electronic supplementary material 4) associated with (i) models based on *K*-theory, (ii) solution with explicit functional forms of *K*_{z}(*x*,*z*), (iii) neglect of longitudinal diffusion, (iv) applicability in low wind speed conditions, (iv) the parameterization of *K*_{z}(*x*,*z*) used in the evaluation and (v) neglecting the veering in wind speed with vertical height.

## 5. Conclusions

For the turbulent dispersion of a pollutant released from a continuous source in the ABL, a generalized analytical model describing the crosswind-integrated concentrations is presented. An analytical scheme is described to solve the resulting two-dimensional steady-state advection–diffusion equation for horizontal wind speed as a generalized function of vertical height above the ground and eddy diffusivity as a function of both downwind distance from the source and vertical height. On the basis of a mathematical analysis, the appropriate functional forms of vertical eddy diffusivity are derived. Particular cases of this model are deduced for (i) uniform wind speed and eddy diffusivity leading to the classical Gaussian plume model and (ii) uniform wind speed and eddy diffusivity as a function of the downwind distance. By taking the particular forms of wind speed and vertical eddy diffusivity, the solutions of the resulting advection–diffusion equation obtained from the present technique are found to be close to those derived using classical techniques. Some of the limitations associated with this work are pointed out.

The model is evaluated with the observations obtained from the Copenhagen diffusion experiment in unstable conditions and the Hanford and Prairie Grass experiments in stable conditions. In evaluating the model, we have used the recently proposed formulation (Gryning *et al.* 2007) for the wind speed in the entire ABL. From the analysis of statistical measures for diffusion experiments in different stability conditions, it is concluded that the present model performs well with the observations.

The present model can also be used for the concentration distribution of a pollutant released from a line source perpendicular to the direction of the mean wind. The model is also useful to study the turbulent dispersion from a steady two-dimensional horizontal source in a wide, open channel for a generalized profile of wind flow and eddy diffusivity as discussed by Nokes *et al.* (1984).

## Appendix A. Evaluation of the forms of *K*_{z}(*x*,*z*)

The commutative condition (2.16) may not be satisfied for all generalized forms of *F*(*x*) corresponding to *K*_{z}(*x*,*z*). To find some of these suitable forms of *K*_{z}(*x*,*z*), we have used the following theorem:

### Theorem (Martin 1967).

*Let F*(*x*) *be the continuous matrix on an interval X. Then the following are equivalent for all x*_{0}*, x*_{1}*and x*_{2}*in X:*
A1
A2
A3

Notice that equation (A3) of the theorem implies that the condition (2.16) is satisfied in a natural way. According to the theorem, it is adequate to verify one of these parts to ensure that the condition (2.16) is satisfied for the matrix *F*(*x*). Based on this theorem, we find the following functional forms of *K*_{z}(*x*,*z*) for matrix *F* satisfying the condition (2.16):
A4
A5
A6
A7
where *f*(*x*) and *g*(*z*) are the generalized functions of *x* and *z*, respectively. In spite of these forms of *K*_{z}(*x*,*z*), there may be other implicit forms of *K*_{z} for matrix *F* satisfying the commutative condition (2.16). Within the analytical framework, a choice of an implicit functional form of *x* and *z* for the eddy diffusivity *K*_{z} is not obvious. However, numerical integration will be needed with an appropriate implicit functional form of *K*_{z} even for evaluating the matrix *F*(*x*). If such an *F*(*x*) satisfies condition (2.16), the solution described in §2*c* can be used for computing the concentration distribution. Most of the parameterizations for *K*_{z}(*x*,*z*) used in atmospheric dispersion modelling are in the forms defined in equations (A4)–(A7) satisfying condition (2.16). Thus, in the present manuscript, we have used these commonly used parameterizations of *K*_{z} to find the solution, in the analytical framework, of the resulting partial differential equation (equation (2.2)) describing the crosswind-integrated concentration.

## Footnotes

- Received July 25, 2009.
- Accepted September 16, 2009.

- © 2009 The Royal Society