## Abstract

In this work, the performance of a unified formal analytical solution for the simulation of atmospheric diffusion problems under stable conditions is evaluated. The eigenquantities required by the formal analytical solution are obtained by solving numerically the associated eigenvalue problem based on a newly developed algorithm capable of being used in high orders and without missing eigenvalues. The performance of the formal analytical solution is evaluated by comparing the converged predicted results against the observed values in the stable runs of the Prairie Grass experiment as well as the simulated results available in the literature. It was found that the developed algorithm was efficient and that the convergence rate depends on the stability condition and the considered parametrizations for wind speed and turbulence. The comparisons among predicted and observed concentrations showed a good agreement and indicate that the considered dispersion formulations are appropriate to simulate dispersion under slightly to moderate atmospheric stable conditions.

## 1. Introduction

Analytical solutions of advective–diffusive transport problems continue to be of interest in many areas of science and engineering, mainly to develop quick and accurate methods for the solution of partial differential equations. A relatively simple model is suitable to carry on sensitivity studies with the objective of evaluating chemical and physical mechanisms as well as to develop physical parametrizations for air quality models involving the closure problem of the turbulence [1].

Efforts have been made to obtain solutions to study the dispersion of a contaminant released from a source in the atmospheric boundary layer (ABL). However, the solutions are mainly available over unstable and neutral conditions [2–5], whereas there is a lack of works considering the solution of the pollutant dispersion equations under stable conditions [3,6,7]. As argued by Salmond & McKendry [8], the strongly layered structure and intermittent characteristics of turbulence in the very stable nocturnal boundary layer present a particular challenge to meteorologists and air pollution modellers alike.

Many processes might occur simultaneously at the stable boundary layer (SBL), indicating the need for a variety of mixing-length scales to represent all the involved physical mechanisms in its description such as buoyancy waves, inertial oscillation, drainage flow and low-level jets [9]. Usually, this layer is non-stationary with no satisfactory proposed theory for the intermittent region yet, and the expressions of profiles of the Reynolds stresses, heat flux and velocity variances are difficult to find [10]. Some aspects of the turbulent SBL structure can be found in the observational study proposed by Mahrt *et al.* [11], Nieuwstadt [12,13] and Holtslag & Nieuwstadt [14].

In atmospheric diffusion modelling, the reliability to represent such phenomena strongly depends on how the turbulence parametrization and micrometeorological parameters are estimated in order to improve the modelling of the physical processes at the ABL [15]. Indeed, there is a compromise between the realistic modelling of the turbulence and the adequate solution of the advection–diffusion equation that represents the pollutant dispersion at the ABL.

Pérez Guerrero *et al.* [4] using the classic integral transform technique (CITT) developed a unified analytical solution of the steady-state atmospheric diffusion equation. The solution was obtained by considering an arbitrary mean wind velocity depending on the vertical coordinate (*z*) and a generalized separable functional form for the eddy diffusivities in terms of the longitudinal (*x*) and vertical coordinates (*z*). The obtained formal analytical solution for a finite domain was given in terms of a series involving eigenvalues, eigenfunctions and norm related to the associated eigenvalue problem of the original diffusive atmospheric equation. This generalized formal analytical solution was extended to semi-infinite and infinite domains by exploring a unified procedure from the directly associated eigenvalue problem. The authors highlighted that the associated eigenvalue problem used in the unified analytical solution improves the convergence behaviour when compared with a simple auxiliary eigenvalue problem.

The solution of the associated eigenvalue problem for situations with arbitrary functional forms for the mean wind velocity and eddy diffusivities was obtained by Pérez Guerrero *et al*. [4] using the generalized integral transform technique (GITT), following the procedure described by Mikhailov & Cotta [16]. With that procedure, the ordinary differential equation of the eigenvalue problem was reduced to an algebraic eigenvalue system. The unified analytical solution developed by Pérez Guerrero *et al*. [4] is valid for all types of atmospheric stability. However, it was applied only for the unstable results of the Prairie Grass and Copenhagen experiments.

The objective of this work is to make an assessment of the performance of the unified formal analytical solution of the atmospheric diffusion equation developed by Pérez Guerrero *et al*. [4] for situations where the transport of pollutants occurs under stable conditions. For this purpose, the stable runs of the classic Prairie Grass experiment [17] are simulated and compared with the observational data of the experiment. A novel algorithm for the solution of the associated eigenvalue problem is developed and extends that suggested by Mikhailov & Ozisik [18] overcoming the risk of missing eigenvalues in higher orders. The versatility/flexibility of the formal unified analytical solution to use eigenquantities obtained regardless of the nature of its algorithms will be shown.

## 2. Problem formulation

Following Pérez Guerrero *et al*. [4], a steady-state concentration distribution *c*≡*c*(*x*,*y*,*z*) of a chemically inert specie is considered, released from a continuous source of strength *Q* located at a specified position (*x*_{s}, *y*_{s}, *z*_{s}). After assuming that the wind flow is aligned with the *x* direction and depends only on the height (*z*), the turbulent diffusion in the direction of the mean wind is neglected compared with the advective transport mechanism and the source is located at *x*_{s}=0, the Eulerian steady-state atmospheric diffusion equation is given by
*u*(*z*), *K*_{y}(*x*,*z*), *K*_{z}(*x*,*z*), respectively, are the wind speed and the transversal and vertical eddy diffusivities; *δ*′ is the Dirac delta generalized function; and *x*, *y* and *z* are the longitudinal, transversal and vertical directions, respectively; [−*y*_{1},*y*_{1}] is the domain in the *y* direction; *z*_{0} is the roughness length; and *z*_{1} is the top of the ABL.

A separable functional structure for the eddy diffusivities is considered as

## 3. Unified formal analytical solution

The solution of equations (2.1)–(2.4) is obtained by applying the systematized CITT procedure given by Ozisik [19], Mikhailov & Ozisik [18] and Cotta [20]. The basic steps can be summarized as follows:

(a) definition of the auxiliary eigenvalue problem,

(b) development of the integral transform pair,

(c) integral transform of the differential equation and

(d) analytical solution for the transformed and original problems.

The formal analytical solution of the diffusion equation (equations (2.1)–(2.4)) for a finite domain was obtained by Pérez Guerrero *et al*. [4] using the CITT. We use the alternative formula given by those authors where the solution is rewritten as the product of two functions, as follows:

Note that the solution was obtained using the CITT and not the solution product procedure. Of course, the solution product procedure could be used for the solution of the atmospheric diffusion equation (equations (2.1)–(2.4)) because the problem is linear, the boundary conditions are homogeneous and the coefficients of the diffusion equation are separable.

The eigenfunctions *Y* _{i}(*y*), eigenvalues *γ*_{i} and norm *N*_{y}(*γ*_{i}) are given by
*et al*. [4] extended the formula of *c*_{y}(*x*,*y*) developed for a finite domain to situations where *y* belongs to an infinite domain, resulting in
*Z*_{j}(*z*) and eigenvalues *η*_{j} are the set of linearly independent solutions related to the following Sturm–Liouville eigenvalue-type problem [4]:
*N*_{z}(*η*_{j}) is associated with the following orthogonality property:
*δ*_{i,j} is the Kronecker delta.

### (a) Solution of the eigenvalue problem in *z*

As observed by Pérez Guerrero *et al*. [4], a closed-form analytical solution for the eigenvalue problem of equations (3.7*a*–*c*) is guaranteed only for special functional forms of *u*(*z*) and *K*(*z*), where the eigenvalue problem can be reduced to the generalized Bessel equation [18,21]. For general functional forms of *u*(*z*) and *K*(*z*), the solution of the eigenvalue problem could be obtained by using the semi-analytic sign-count method [22] or from algorithms for automatic computation of eigenvalue and eigenfunctions of Sturm–Liouville problems [23,24]. The procedure using the GITT given by Cotta [20], and Mikhailov & Cotta [16], is another alternative to solve the eigenvalue problem. In this case, the eigenvalue differential equation is transformed in an algebraic eigenvalue system, which is then solved numerically or analytically such as described in Cotta [20]. Finding the eigenvalues of a Sturm–Liouville problem can be a computationally challenging task, especially when a large set of eigenvalues is computed, or just when particularly large eigenvalues are sought. This is a consequence of the highly oscillatory behaviour of the solutions corresponding to high eigenvalues, which forces a naive integrator to take increasingly smaller steps [25]. In this context, such as reported by Mikhailov & Ozisik [18], Cotta [20] and recently by Pérez Guerrero *et al*. [26] solving the advection–dispersion transport equation in layered media, a difficulty in the numerical solution of an eigenvalue problem is related to the risk of missing some of the eigenvalues.

Mikhailov & Ozisik [18] used an algorithm with standard Runge–Kutta and Newton's iterative method to calculate numerically the eigenvalues, eigenfunctions and norm of the Sturm–Liouville system. About that algorithm, those authors observed that: ‘the method is convenient for computing the lower eigenvalues, but difficulties may arise in calculating the higher order ones because of the rapid oscillation character of the higher eigenfunctions. When using such algorithms, extreme care must be exercised to ensure that no eigenvalues are missed’.

Motivated to overcome the limitations in observations given by Mikhailov & Ozisik [18] in their algorithm, we opted for a numerical solution of the differential eigenvalue problem with the intention of highlighting the versatility/flexibility of the formal unified analytical solution in use of eigenquantities obtained regardless of the nature of its algorithms. The interesting and novel feature in our numeric calculation of the eigenquantities is the modification of the algorithm given by Mikhailov & Ozisik [18]. Differently from that of those authors, our approach uses an integrator of ordinary differential equations for initial value problems capable of treating stiff situations instead of the standard Runge–Kutta method. In actuality, general-purpose integrators of stiff ordinary differential equations with user-prescribed accuracy exist. Among them can be mentioned the NDSolve from Mathematica software [27], which is a general purpose subroutine capable of solving stiff ordinary differential equations. The algorithms of the literature dedicated to solve the Sturm–Liouville problems do not use directly general-purpose integrators.

On the other hand, another crucial aspect for the numeric solution of the eigenvalue problem is the initial guess. An inappropriate initial guess will give a duplication of eigenvalues calculated or the risk of missing any one of them. The initial guess can be developed from the asymptotic behaviour of the Sturm–Liouville equation. The asymptotic formula for the *i*th eigenvalue is given by Mikhailov & Ozisik [18] as
*ε*_{i} represents a small number associated with the *i*th eigenvalue *η*_{i}.

Mikhailov & Ozisik [18] suggested, as an initial guess, equation (3.9) with *ε*_{i}=0. However, from several tests it was verified that this suggestion was not appropriate because of the duplication and missing of calculated eigenvalues. We developed an expression for the initial guess. Thus, considering *ε*_{i+1}≈*ε*_{i}, the (*i*+1)th eigenvalue can be estimated from equation (3.9), resulting in
*b*,*c*), *η*_{0}=0 is the first eigenvalue with eigenfunction *Z*_{0}(*z*)=1. In this way, the original algorithm given by Mikhailov & Ozisik [18] was modified considering a general purpose subroutine capable of solving stiff ordinary differential equations and using equation (3.10) as an initial guess for the (*i*+1)th eigenvalue. In the algorithm, the number of iterations desired for the determination of the eigenvalue is also defined and the relative and absolute error for each iteration are calculated.

## 4. Results and discussion

The performance of the formal unified analytical solution given in §3 is evaluated by simulating only the stable conditions of the classic Prairie Grass experiment carried out in O'Neill, NE, USA, in July and August, 1956 [17]. In the experiment, SO_{2} was released continuously (10 min) from a near-surface source (almost all runs at 0.46 m and in some of them at 1.5 m), under a variety of meteorological conditions. The site was a flat and homogeneous prairie covered with short grass. The sampling of the released material was done at a height of 1.5 m on five concentric arcs located 50, 100, 200, 400 and 800 m downwind from the source.

Tables 1 and 2 summarize the micrometeorological parameters, emission data and the observed crosswind-integrated concentrations normalized by the source strength. The very stable cases (*z*_{1}/*L*>10) were not considered, based on the same argument given by Kumar & Sharan [3], i.e. that the turbulence becomes very weak, sporadic, intermittent and no longer continuous in time and space [14]. For comparative purposes, we use the same digital version of the dataset used by Kumar & Sharan [3] as available at http://www2.dmu.dk/atmosphericenvironment/Docs/PrairieGrass.xls, that gives the normalized crosswind-integrated concentrations and meteorological variables with a surface roughness length *z*_{0}=0.006 m.

The simulations were done considering the following two different configurations, CONFIG1 and CONFIG2. On the first one, the wind velocity and momentum eddy diffusivity profiles are such as proposed by Ulke [28]. On the second one, the velocity profile is given by Gryning *et al*. [29] and the eddy diffusivity by Degrazia & Moraes [30], as adopted by Kumar & Sharan [3]. It is important to note that Ulke's parametrizations consider the current understanding of the ABL in a simple and continuous formulation through the different atmospheric conditions, where the stability-dependent function provides a smooth variation between the different regimes in the ABL, in accordance with Holstlag and Nieuwstadt's proposal based on the relevant key variables. Besides this, it includes the effects of mechanical and buoyancy-induced/suppressed turbulence in a consistent way, in the wind and the eddy diffusivity functional forms [4].

The quantitative evaluation was done using the analysis of the convergence behaviour of the results and scatter diagrams for the observed and predicted ground-level crosswind-integrated concentrations. The ability of the modelling approach to estimate the whole range of measured concentrations was studied using residual plots. The statistical indices proposed by Hanna *et al*. [31] were used to analyse the agreement between the observed (*O*) and predicted (*P*) values. The indices fractional bias (FB), fractional variance (FS), normalized mean square error (NMSE), correlation coefficient (*R*) and fraction within a factor of 2 (FAC2) are defined as follows:
*σ*_{O} and *σ*_{P} are the standard deviations of observed and predicted quantities, respectively, and the over-bar indicates an average.

In the following sections, the performance of the formal analytical solution given in §3 is discussed from the analysis of the convergence of equation (3.3) and by evaluating the simulated crosswind-integrated concentrations and compared against those obtained from the observed values in the experiment and simulated results from the open literature.

### (a) Calculation of the eigenquantities

The algorithm for the calculation of the eigenquantities was implemented in the Mathematica software [27], where 13 iterations were specified. The final maximum relative and absolute errors for the determination of the eigenvalues were typically in the order of 10^{−8} and 10^{−9}. Tables 3 and 4 show the typical number of iterations for the determination of the converged eigenvalues for run 68 of the Prairie Grass experiment using both CONFIG1 and CONFIG2. The results indicate a fast convergence in the calculation of the eigenvalues, where only four iterations are required for the first eigenvalue, decreasing the number of iterations for the second and tenth eigenvalues. Intermediate or high eigenvalues require only one iteration to obtain converged eigenvalues. This fast convergence is the effect of the initial guess used (equation (3.10)). Although not shown here, other simulations were also done to evaluate the performance of the algorithm for the calculation of the eigenquantities, including very high eigenvalues (*i*=500).

### (b) Convergence analysis of the simulated results

Experiments 24 and 68 were chosen to illustrate the convergence of the solution because they represent two extreme situations given by a slightly stable condition with *z*_{1}/*L*=0.61 (run 24) and a moderately stable condition with *z*_{1}/*L*=8.50 (run 68).

Figure 1*a*,*b* shows the convergence behaviour of the simulated normalized crosswind-integrated concentration *c*_{z}(*x*,1.5)/*Q* (s m^{−2}) for a different number of summed terms (*N*) using the formal analytical solution (equation (3.3)). Typically, for all the runs of the simulated experiment the number of summed terms in the series of equation (3.3) increases as the longitudinal position is closer to the source. On the other hand, from figure 1*a*,*b* it is also noted that the convergence rate is influenced by the stability condition and by the used parametrizations. For slightly stable conditions (run 24), shown in figure 1*a*, the convergence of the results obtained using parametrizations defined by CONFIG1 is slightly slower than that obtained with CONFIG2. Thus, calculation of converged results with a high accuracy of six significant figures for the slightly stable conditions of run 24 required 30–90 summed terms at the different sampling positions of the experiment using parametrizations of CONFIG1, while 20–80 summed terms were required for the same positions when CONFIG2 is used. However, for a moderate stable condition (run 68), shown in figure 1*b*, the convergence rate of the results with the parametrizations defined in CONFIG1 is much slower than that of CONFIG2. Thus, run 68 required 80–260 summed terms to reach the converged results with six significant figures when CONFIG2 is used, while 120–400 summed terms were required to obtain the same significant figures in the converged results when CONFIG1 is adopted. For engineering purposes, results with an accuracy of two significant figures can be obtained with less than *N*=120 summed terms.

The different numbers of terms required to achieve convergence for CONFIG1 and CONFIG2 can be explained from the influence of the diffusive and the advective terms in the atmospheric diffusion equation. According to Cotta [20] and Pérez Guerrero *et al*. [32], a faster convergence corresponds to a smaller Peclet (Pe) number. This can be verified defining *Pe*=1025 (CONFIG1) and 806 (CONFIG2); and for run 68: *Pe*=22267 (CONFIG1) and 11189 (CONFIG2). Those values of the Pe number explain the faster convergence of CONFIG2 in relation to CONFIG1.

### (c) Analysis of the simulated results

Figure 2 presents the converged predicted values of the normalized crosswind-integrated concentrations for the CONFIG1 and CONFIG2 cases against Prairie Grass experimental data. The comparison shows a good agreement between the modelled values and the experimental ones, further indicating that the simulated results using CONFIG 1 and CONFIG 2 are relatively close and have a similar trend. The simulated results for both cases (CONFIG1 and CONFIG2) tend, in general, to overestimate the experimental data. Thus, at positions *x*=200, 400 and 800 m all the results were overestimated in relation to the observed values, whereas for *x*=50 and 100 m some results overestimate and others underestimate the experimental values. It must be mentioned that Kumar & Sharan [3] obtained similar results. These effects of overestimation and underestimation of the experimental values can be appreciated in the scatter diagram of figure 2, where the higher proportion of overestimated values occurs for the smaller normalized concentrations. However, almost the majority of observed crosswind-integrated concentrations are predicted within a factor of 2.

Figure 3 illustrates the variation of the ratio of predicted to observed normalized crosswind-integrated concentrations, with the observed ones using both options for the formulation of the tracer dispersion. Almost all the ratio values are near the value of 1, which indicates a very good agreement. The overprediction of the predicted concentration values for positions farther from the source may be related to the limitations of the used parametrizations in representing the real mechanisms in the atmospheric diffusion process under stable conditions. In addition, the overestimation might also be related to the dry deposition of the tracer that was not accounted for in the mathematical model.

Figures 4 and 5 present the converged values of the normalized crosswind concentrations for the CONFIG1 and CONFIG2 cases against results from the Operationelle Meteorologiske Luftkvalitetsmodeller (OML) model proposed by Olesen *et al.* [33] available at http://www2.dmu.dk/atmosphericenvironment/Docs/PrairieGrass.xls. The comparison shows a good agreement between the modelling results, confirming the performance of the developed model in both configurations.

The statistical measures obtained in the evaluation of the performance of the formal analytical solution for cases with CONFIG1 and CONFIG2 are presented in table 5. The comparison was done with other numerical results proposed by Moreira *et al.* [6], Kumar & Sharan [3] and the OML model. Moreira *et al.* [6] developed a Eulerian model and a Lagrangian particle model to simulate the dispersion of a contaminant released from a low-level source in the SBL using two different eddy diffusion parametrizations proposed by Hanna [34] and Degrazia *et al*. [35]. The friction velocity and Monin–Obukhov length were based on the paper of van Ulden [36], and the formula proposed by Zilitinkevich [37] was used to calculate the ABL height. Kumar & Sharan [3] proposed an analytical scheme to solve the resulting two-dimensional steady-state advection–diffusion equation. For the stable runs of the Prairie Grass experiment, the simulated concentrations were obtained using the eddy diffusion parametrizations proposed by Degrazia & Moraes [30] and the extended profile of wind speed into the entire ABL given by Gryning *et al*. [29].

For both sets of parametrizations used in the proposed model, the statistical indices are similar and with values close to the ideal ones, with a slightly overall better performance for CONFIG1. Therefore, the simulation of the considered stable Prairie Grass experimental runs with the parametrizations defined by CONFIG1 and CONFIG2 is appropriate to simulate slightly and moderately stable situations. The indices obtained with CONFIG1 and CONFIG2 are better than those found by Kumar & Sharan [3] and Moreira *et al.* [6]. Considering that CONFIG2 was also used by Kumar & Sharan [3], the difference in the results can be explained by the use of an associated eigenvalue problem for the determination of the eigenquantities (eigenfunction, eigenvalue and norm) in our formal analytical solution while Kumar & Sharan [3] used a simple eigenvalue problem. The discrepancy of our results with the Moreira *et al.* [6] could be related to the different parametrization sets used. Moreira *et al.* [6] considered the parametrization suggested by Degrazia *et al*. [35], which is different from that defined in CONFIG1 and CONFIG2.

The FB indexes indicate that observed concentrations are overestimated in the mean by all modelling approaches with the greatest bias for that of Moreira *et al.* [6] and the smallest for the present model with CONFIG1. The FSs show that the models underpredict the scattering of the measured concentrations, with the exception of the Moreira *et al.* [6] model. The latter shows also the greatest variability, with an NMSE an order of magnitude higher than the other models. The correlation coefficients show a high degree of association between the observed and predicted concentrations, with the greatest values for Kumar & Sharan [3], followed by those of the present model.

## 5. Conclusion

The performance of the unified formal analytical solution of the atmospheric diffusion equation developed by Pérez Guerrero *et al*. [4] for situations where the transport of pollutants occurs under stable conditions was evaluated with the Prairie Grass data for experimental runs with *z*_{1}/*L*<10 and results from other models available in the open literature.

The newly developed algorithm for the numerical solution of the eigenvalue problem in the calculation of the eigenquantities used in the formal analytical solution was efficient, capable of being used in higher order and without missing eigenvalues.

The comparison with modelling approaches available in the literature highlights the potentialities of the proposed unified solution, on the one hand indicating the benefits in terms of the reduction of the global measures of error, and, on the other hand, as a procedure to test the capability of different dispersion parametrizations. In this way, the present model has proved its suitability for analysing the performance of proposed formulations for dispersion in stable conditions.

The knowledge of dispersion under stable conditions is still a challenging issue, particularly linked to the availability of experimental data to evaluate the proposed formulations. In the SBL, which is rarely stationary, the turbulence is often intermittent and patchy and strongly depends on the flow conditions. The measurements under such conditions are hard to carry out. Most of the applied parametrizations are based on assumptions and hypotheses that are difficult to find in a real SBL.

## Funding statement

The authors acknowledge the support of the Brazilian National Research Council (CNPq) and Carlos Chagas Filho Foundation of Rio de Janeiro State (FAPERJ), and the partial support of the projects UBACyT: 20020100101013 financed by the University of Buenos Aires and PICT08-1739 from the National Agency for the Promotion of Science and Technology (ANPCyT), Argentina.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.