Analytical Model for Air Pollution in the Atmospheric Boundary Layer

The worldwide concern due to the increasing frequency of ecological disasters on our planet urged the scientific community to take action. One of the possible measures are analytical descriptions of pollution related phenomena, or simulations by effective and operational models with quantitative predictive power. The atmosphere is considered the principal vehicle by which pollutant materials are dispersed in the environment, that are either released from the productive and private sector or eventually in accidental events and thus may result in contamination of plants, animals and humans. Therefore, the evaluation of airborne material transport in the Atmospheric Boundary Layer (ABL) is one of the requirements for maintenance, protection and recovery of the ecological system. In order to analyse the consequences of pollutant discharge, atmospheric dispersion models are of need, which have to be tuned using specific meteorological parameters and conditions for the considered region. Moreover, they shall be subject to the local orography and shall supply with realistic information concerning environmental consequences and further help reduce impact from potential accidents such as fire events and others. Moreover, case studies by model simulations may be used to establish limits for the escape of pollutant material from specific sites into the atmosphere.


Introduction
The worldwide concern due to the increasing frequency of ecological disasters on our planet urged the scientific community to take action. One of the possible measures are analytical descriptions of pollution related phenomena, or simulations by effective and operational models with quantitative predictive power. The atmosphere is considered the principal vehicle by which pollutant materials are dispersed in the environment, that are either released from the productive and private sector or eventually in accidental events and thus may result in contamination of plants, animals and humans. Therefore, the evaluation of airborne material transport in the Atmospheric Boundary Layer (ABL) is one of the requirements for maintenance, protection and recovery of the ecological system. In order to analyse the consequences of pollutant discharge, atmospheric dispersion models are of need, which have to be tuned using specific meteorological parameters and conditions for the considered region. Moreover, they shall be subject to the local orography and shall supply with realistic information concerning environmental consequences and further help reduce impact from potential accidents such as fire events and others. Moreover, case studies by model simulations may be used to establish limits for the escape of pollutant material from specific sites into the atmosphere.
The theoretical approach to the problem essentially assumes two basic forms. In the Eulerian approach, diffusion is considered, at a fixed point in space, proportional to the local gradient of the concentration of the diffused material. Consequently, it considers the motion of fluid within a spatially fixed system of reference. These kind of approaches are based on the resolution, on a fixed space-time grid, of the equation of mass conservation of the pollutant chemical species. Lagrangian models are the second approach and they differ from Eulerian ones in adopting a system of reference that follows atmospheric motions. Initially, the term Lagrangian was used only to refer to moving box models that followed the mean wind trajectory. Currently, this class includes all models that decompose the pollutant cloud into discrete "elements", such as segments, puffs or pseudo-particles (computer particles). In Lagrangian particle models pollutant dispersion is simulated through the motion of computer 2 www.intechopen.com particles whose trajectories allow the calculation of the concentration field of the emitted substance. The underlying hypothesis is that the combination of the trajectories of such particles simulates the paths of the air particles situated, at the initial instant, in the same position. The motion of the particles can be reproduced both in a deterministic way and in a stochastic way.
In the present discussion we adopt an Eulerian approach, more specifically to the K model, where the flow of a given field is assumed to be proportional to the gradient of a mean variable. K-theory has its own limits, but its simplicity has led to a widespread use as the mathematical basis for simulating pollution dispersion. Most of of Eulerian models are based on the numerical resolution of the equation of mass conservation of the pollutant chemical species. Such models are most suited to confronting complex problems, for example, the dispersion of pollutants over complex terrain or the diffusion of non-inert pollutants.
In the recent literature progressive and continuous effort on analytical solution of the advection-diffusion equation (ADE) are reported. In fact analytical solutions of equations are of fundamental importance in understanding and describing the physical phenomena. Analytical solutions explicitly take into account all the parameters of a problem, so that their influence can be reliably investigated and it is easy to obtain the asymptotic behaviour of the solution, which is usually more tedious and time consuming when generated by numerical calculations. Moreover, like the Gaussian solution, which was the first solution of ADE with the wind and the K diffusion coefficients assumed constant in space and time, they open pathways for constructing more sophisticated operative analytical models. Gaussian models, so named because they are based on the Gaussian solution, are forced to represent real situations by means of empirical parameters, referred to as "sigmas". In general Gaussian models are fast, simple, do not require complex meteorological input, and describe the diffusive transport in an Eulerian framework exploring the use of the Eulerian nature of measurements. For these reasons they are widely employed for regulatory applications by environmental agencies all over the world although their well known intrinsic limits. For the further discussion we restrict our considerations to scenarios in micro-meteorology were the Eulerian approach has been validated against experimental data and has been found useful for pollution dispersion problems in the ABL. A significant number of works regarding ADE analytical solution (mostly two-dimensional solutions) is available in the literature. By analytical we mean that starting from the advection-diffusion equation, which we consider adequate for the problem, the solution is derived without introducing approximations that might simplify the solution procedure. A selection of references considered relevant by the authors are Rounds (1955), Smith (1957), Scriven & Fisher (1975, Demuth (1978), van Ulden (1978, Nieuwstadt & Haan (1981), Tagliazucca et al. (1985), Tirabassi (1989), Tirabassi & Rizza (1994), Sharan et al. (1996), Lin & Hildemann (1997), Tirabassi (2003). It is noteworthy that solutions from citations above are valid for specialized problems. Only ground level sources or infinite height of the ABL or specific wind and vertical eddy diffusivity profiles are admitted. A more general approach that may be found in the literature is based on the Advection Diffusion Multilayer Method (ADMM), which solves the two-dimensional ADE with variable wind profile and eddy diffusivity coefficient (Moreira et al. (2006a)). The main idea here relies on the discretisation of the Atmospheric Boundary Layer (ABL) in a multilayer domain, assuming in each layer that the eddy diffusivity and wind profile take averaged values. The resulting advection-diffusion equation in each layer is then solved by the Laplace Transform technique. For more details about this methodology see reference Moreira et al. (2006a).
A further mile-stone in the direction of the present approach is given in references Costa et al. (2006; by the Generalized Integral Advection Diffusion Multilayer Technique (GIADMT), that presents a general solution for the time-dependent three-dimensional ADE again assuming the stepwise approximation for the eddy diffusivity coefficient and wind profile and proceeding further in a similar way as the multi-layer approach (Moreira et al. (2006a)), however in a semi-analytical fashion. In this chapter we report on a generalisation of the afore mentioned models resulting in a general space-time (3 ⊕ 1) solution for this problem, given in closed form.
To this end, we start from the three-dimensional and time dependent ADE for the ABL as the underlying physical model. In the subsequent sections we show how the model is solved analytically using integral transform and a spectral theory based methods which then may provide short, intermediate and long term (normalized) concentrations and permit to assess the probability of occurrence of high contamination level case studies of accidental scenarios. The solution of the 3 ⊕ 1 space-time solution is compared to other dispersion process approaches and validated against experiments. Comparison between different approaches may help to pin down computational errors and finally allows for model validation. In this line we show with the present discussion, that our analytical approach does not only yield an acceptable solution for the time dependent three dimensional ADE, but further predicts tracer concentrations closer to observed values compared to other approaches from the literature. This quality is also manifest in better statistical coefficients including model validation.
Differently than in most of the approaches known from the literature, where solutions are typically valid for scenarios with strong restrictions with respect to their specific wind and vertical eddy diffusivity profiles, the discussed method is general in the sense that once the model parameter functions are determined or known from other sources, the hierarchical procedure to obtain a closed form solution is unique. More specifically, the general analytical solution for the advection-diffusion problem works for eddy diffusivity and wind profiles, that are arbitrary functions with a continuous dependence on the vertical and longitudinal spatial variables, as well as time.
A few technical details are given in the following, that characterise the principal steps of the procedure. In a formal step without physical equivalence we expand the contaminant concentration in a series in terms of a set of orthogonal eigenfunctions. These eigenfunctions are the solution of a simpler but similar problem to the existing one and are detailed in the next section. Replacing this expansion in the time-dependent, three-dimensional ADE in Cartesian geometry, we project out orthogonal components of the equation, thus inflating the tree-dimensional advection-diffusion equation into an equation system. Note, that the projection integrates out one spatial dimension, so that the resulting problem reduces the problem effectively to 2 ⊕ 1 space-time dimensions. Such a problem was already solved by the Laplace transform technique as shown in Moreira et al. (2009a), Moreira et al. (2009b), Buske et al. (2010) and references therein. In the next sections we present the prescription for constructing the solution for the general problem, but for convenience and without restricting generality we resort to simplified models. Further, we present numerical simulations and indicate future perspectives of this kind of modelling.

41
Analytical Model for Air Pollution in the Atmospheric Boundary Layer www.intechopen.com

The advection-diffusion equation
The advection-diffusion equation of air pollution in the atmosphere is essentially a statement of conservation of the suspended material and it can be written as (Blackadar (1997)) wherec denotes the mean concentration of a passive contaminant (in units of g m 3 ), D Dt is the substantial derivative, ∂ t is a shorthand for the partial time derivative andŪ =(ū,v,w) T is the mean wind (in units of m s ) with Cartesian components in the directions x, y and z, respectively, ∇ =( ∂ x , ∂ y , ∂ z ) T is the usual Nabla symbol and S is the source term. The terms U ′ c ′ represent the turbulent flux of contaminants (in units of g sm 2 ), with its longitudinal, crosswind and vertical components.
Observe that equation (1) has four unknown variables (the concentration and turbulent fluxes) which lead us to the well known turbulence closure problem. One of the most widely used closures for equation (1), is based on the gradient transport hypothesis (also called K-theory) which, in analogy with Fick's law of molecular diffusion, assumes that turbulence causes a net movement of material following the negative gradient of material concentration at a rate which is proportional to the magnitude of the gradient (Seinfeld & Pandis (1998)).
Here, the eddy diffusivity matrix K = diag(K x , K y , K z ) (in units of m 2 s ) is a diagonal matrix with the Cartesian components in the x, y and z directions, respectively. In the first order closure all the information on the turbulence complexity is contained in the eddy diffusivity.
The simplicity of the K-theory of turbulent diffusion has led to the widespread use of this theory as a mathematical basis for simulating pollutant dispersion (open country, urban, photochemical pollution, etc.), but K-closure has its known limits. In contrast to molecular diffusion, turbulent diffusion is scale-dependent. This means that the rate of diffusion of a cloud of material generally depends on the cloud dimension and the intensity of turbulence. As the cloud grows, larger eddies are incorporated in the expansion process, so that a progressively larger fraction of turbulent kinetic energy is available for the cloud expansion.
Equation (3) is considered valid in the domain (x, y, z) ∈ Γ bounded by 0 < x < L x ,0< y < L y and 0 < z < h and subject to the following boundary and initial conditions, Instead of specifying the source term as an inhomogeneity of the partial differential equation, we consider a point source located at an edge of the domain, so that the source position r S = (0, y 0 , H S ) is located at the boundary of the domain r S ∈ δΓ. Note, that in cases where the source is located in the domain, one still may divide the whole domain in sub-domains, where the source lies on the boundary of the sub-domains which can be solved for each sub-domain separately. Moreover, a set of different sources may be implemented as a superposition of independent problems. Since the source term location is on the boundary, in the domain this term is zero everywhere (S(r) ≡ 0 for r ∈ Γ\δΓ), so that the source influence may be cast in form of a condition, where we assume that our coordinate system is oriented such that the x-axis is aligned with the mean wind direction. Since the flow crosses the plane perpendicular to the propagation (here the y-z-plane) the source condition reads where Q is the emission rate (in units of

A closed form solution
In this section we first introduce the general formalism to solve a general problem and subsequently reduce the problem to a more specific one, that is solved and compared to experimental findings.

The general procedure
In order to solve the problem (3) we reduce the dimensionality by one and thus cast the problem into a form already solved in reference Moreira et al. (2009a). To this end we apply the integral transform technique in the y variable, and expand the pollutant concentration as where R =( R 1 , R 2 ,...) T and Y =( Y 1 , Y 2 ,...) T is a vector in the space of orthogonal eigenfunctions, given by Y m (y)=cos(λ m y) with eigenvalues λ m = m π L y for m = 0, 1, 2, . . .. For convenience we introduce some shorthand notations, ∇ 2 =( ∂ x ,0,∂ y ) T and∂ y = (0, ∂ y ,0) T , so that equation (3) reads now, Upon application of the integral operator here F is an arbitrary function and ∧ signifies the dyadic product operator, and making use of orthogonality renders equation (8) a matrix equation. The appearing integral terms are Here, B 0 = L y 2 I, where I is the identity, the elements (Z) mn = 2 1−n 2 /m 2 δ 1,j with δ i,j the Kronecker symbol and j =( m + n)mod2 is the remainder of an integer division (i.e. is one for m + n odd and zero else). Note, that the integrals Ω i and T i depend on the specific form of the eddy diffusivity K. The integrals (10) are general, but for practical purposes and for application to a case study we truncate the eigenfunction space and consider M components in R and Y only, though continue using the general nomenclature that remains valid. The obtained matrix equation determines now together with initial and boundary condition uniquely the components R i for i = 1, . . . , M following the procedure introduced in reference Moreira et al. (2009a).

A specific case for application
In order to discuss a specific case we recall the convention already adopted above, that considered the average wind velocity aligned with the x-axis. Since we consider length scales L x , L Y , h typical for micro-meteorological scenarios we simplifyŪ =(ū,0,0) T . By comparison of physically meaningful cases, one finds for the operator norm ||∂ x K x ∂ x || << |ū|, which can be understood intuitively because eddy diffusion is observable predominantly perpendicular to the mean wind propagation. As a consequence we neglect the terms with K x and ∂ x K x .
The principal aspect of interest in pollution dispersion is the vertical concentration profile, that responds strongly to the atmospheric boundary layer stratification, so that the simplified eddy diffusivity K → K 1 = diag(0, K y , K z ) depends in leading order approximation only on the vertical coordinate K 1 = K 1 (z). For this specific case the integrals Ω i reduce to T 1 → 0 , where Λ = diag(λ 2 1 , λ 2 2 ,...). The simplified equation system to be solved is then, which is equivalent to the problem by virtue of B being a diagonal matrix.
The specific form of the eddy diffusivity determines now whether the problem is a linear or non-linear one. In the linear case the K is assumed to be independent ofc, whereas in more realistic cases, even if stationary, K may depend on the contaminant concentration and thus renders the problem non-linear. However, until now no specific law is known that links the eddy diffusivity to the concentration so that we hide this dependence using a phenomenologically motivated expression for K which leaves us with a partial differential equation system in linear form, although the original phenomenon is non-linear. In the example below we demonstrate the closed form procedure for a problem with explicit time dependence, which is novel in the literature.
The solution is generated making use of the decomposition method (Adomian (1984;1988;) which was originally proposed to solve non-linear partial differential equations, followed by the Laplace transform that renders the problem a pseudo-stationary one. Further we rewrite the vertical diffusivity as a time average termK z (z) plus a term representing the variations κ z (z, t) around the average for the time interval of the measurement K z (x, z, t)= K z (z)+κ z (z, t) and use the asymptotic form of K y , which is then explored to set-up the structure of the equation that defines the recursive decomposition scheme.
The function R = ∑ j R j = 1 T R (c) is now decomposed into contributions to be determined by recursion. For convenience we introduced the one-vector 1 =( 1, 1, . . .) T and inflate the vector R to a vector with each element being itself a vector R j . Upon inserting the expansion in equation (16) one may regroup terms that obey the recursive equations and starts with the time averaged solution for K z .
The extension to the closed form recursion is then given by From the construction of the recursion equation system it is evident that other schemes are possible. The specific choice made here allows us to solve the recursion initialisation using the procedure described in references Moreira et al. (2006a;2009a), where a stationary K was assumed. For this reason the time dependence enters as a known source term from the first recursion step on.

Recursion initialisation
The boundary conditions are now used to uniquely determine the solution. In our scheme the initialisation solution that contains R 0 satisfies the boundary conditions (equations (4)- (6)) while the remaining equations satisfy homogeneous boundary conditions. Once the set of problems (18) is solved by the GILTT method, the solution of problem (3) is well determined.
It is important to consider that we may control the accuracy of the results by a proper choice of the number of terms in the solution series.
In references Moreira et al. (2009a;, a two dimensional problem with advection in the x direction in stationary regime was solved which has the same formal structure than (18) except for the time dependence. Upon rendering the recursion scheme in a pseudo-stationary form problem and thus matching the recursive structure of Moreira et al. (2009a;, we apply the Laplace Transform in the t variable, (t → r) obtaining the following pseudo-steady-state problem.
The x and z dependence may be separated using the same reasoning as already introduced in (7). To this end we pose the solution of problem (19) in the form: where C =( ζ 1 (z), ζ 2 (z),...) T are a set of orthogonal eigenfunctions, given by ζ i (z)= cos(γ l z), and γ i = iπ/h (for i = 0, 1, 2, . . .) are the set of eigenvalues.
Replacing equation (20) in equation (19) and using the afore introduced projector (9) now for the z dependent degrees of freedom h 0 dzC[F]= h 0 F T ∧ C dz yields a first order differential equation system.
where P 0 = P 0 (x, r) and B −1 1 B 2 come from the diagonalisation of the problem. The entries of matrices B 1 and B 2 are A similar procedure leads to the source condition for (21). whereP(s, r) denotes the Laplace Transform of P(x, r). Here X (−1) is the (inverse) matrix of the eigenvectors of matrix B −1 1 B 2 with diagonal eigenvalue matrix D and the entries of matrix (sI + D) ii = s + d i . After performing the Laplace transform inversion of equation (23), we come out with P(x, r)=XG(x, r)X −1 Ω , where G(x, r) is the diagonal matrix with components (G) ii = e −d i x . Further the still unknown arbitrary constant matrix is given by Ω = X −1 P(0, r).
The analytical time dependence for the recursion initialisation (19) is obtained upon applying the inverse Laplace transform definition.
To overcome the drawback of evaluating the line integral appearing in the above solution, we perform the calculation of this integral by the Gaussian quadrature scheme, which is exact if the integrand is a polynomial of degree 2M − 1inthe 1 r variable.
where a and p are respectively vectors with the weights and roots of the Gaussian quadrature scheme (Stroud & Secrest (1966)), and the argument (x, z, p t ) signifies the k-th component of p in the k-th row of pR 0 . Note, k is a component from contraction with a. Concerning the issue of possible Laplace inversion schemes, it is worth mentioning that this approach is exact, however we are aware of the existence of other methods in the literature to invert the Laplace transformed functions , ), but we restrict our attention in the considered problem to the Gaussian quadrature scheme, which is sufficient in the present case. Finally, the solution of the remaining equations of the recursive system (18) are attained in an analogue fashion expanding the source term in series and solving the resulting linear first order differential matrix equation with known and integrable source by the Laplace transform technique as shown above.

Reduction to solutions of simpler models
To establish the solution of the stationary, three and two-dimensional, ADE under Fickian flow regime, with and without longitudinal diffusion, one only need to take the limit t → ∞ of the previously derived solutions, which is equivalent to take r → 0 (Buske et al. (2007a)). Further models that make use of simplifications are not discussed in this chapter because they can be determined with the present developments and applying them to the models of references

Solution validation against experimental data
The solution procedure discussed in the previous section was coded in a computer program and will be available in future as a program library add-on. In order to illustrate the suitability of the discussed formulation to simulate contaminant dispersion in the atmospheric boundary layer, we evaluate the performance of the new solution against experimental ground-level concentrations.

Time dependent vertical eddy diffusivity
As already mentioned in the derivation of the solution, one has to make the specific choice for the turbulence parametrisation. From the physical point of view the turbulence parametrisation is an approximation to nature in the sense that we use a mathematical model as an approximated (or phenomenological) relation that can be used as a surrogate for the natural true unknown term, which might enter into the equation as a non-linear contribution. The reliability of each model strongly depends on the way turbulent parameters are calculated and related to the current understanding of the ABL (Degrazia (2005)).
The present parametrisation is based on the Taylor statistical diffusion theory and a kinetic energy spectral model. This methodology, derived for convective and moderately unstable conditions, provides eddy diffusivities described in terms of the characteristic velocity and length scales of energy-containing eddies. The time dependent eddy diffusivity has been derived by Degrazia et al. (2002) and can be expressed as the following formula.
Here α = x, y, z are the Cartesian components, i = u, v, w are the longitudinal and transverse wind directions, c i = α i (0.5±0.05)(2πκ) − 2 3 , α i = 1 and 4/3 for u, v and w components, respectively (Champagne et al. (1977)), κ = 0.4 is the von Karman constant, ( f * m ) i is the normalized frequency of the spectral peak, h is the top of the convective boundary layer height, w * is the convective velocity scale, ψ is the non-dimensional molecular dissipation rate function and is the non-dimensional travel time, where u is the horizontal mean wind speed.
For horizontal homogeneity the convective boundary layer evolution is driven mainly by the vertical transport of heat. As a consequence of this buoyancy driven ABL, the vertical dispersion process of scalars is dominant when compared to the horizontal ones. Therefore, the present analysis will focus on the vertical time dependent eddy diffusivity only. This vertical eddy diffusivity can be derived from equation (27) by assuming that c w = 0.36 and for the vertical component (Degrazia et al. (2000)). Furthermore, considering the horizontal plane, we can idealize the turbulent structure as a homogeneous one with the length scale of energy-containing eddies being proportional to the convective boundary layer height h.  Gryning & Lyck (1984).
Moreover, for the lateral eddy diffusivity we used the asymptotic behaviour of equation (27) for large diffusion travel times with c v = 0.36 and ( f * m ) v = 0.66 z h (Degrazia et al. (1997)). The dissipation function ψ according to Hφjstrup (1982) and Degrazia et al. (1998) has the form where L is the Obukhov length in the surface layer.

Numerical results
The measurements of the contaminant dispersion in the atmospheric boundary layer consist typically from a sequence of samples over a time period. The experiment used to validate the previously introduced solution was carried out in the northern part of Copenhagen and is described in detail by Gryning & Lyck (1984). Several runs of the experiment with changing meteorological conditions were considered as reference in order to simulate time dependent contaminant dispersion in the boundary layer and to evaluate the performance of the discussed solutions against the experimental centreline concentrations.
The essential data of the experiment are reported in the following. This experiment consisted of a tracer released without buoyancy from a tower at a height of 115m, and a collection of tracer sampling units were located at the ground-level positions up to the maximum of three crosswind arcs. The sampling unit distances varied between two to six kilometers from the point of release. The site was mainly residential with a roughness length of the 0.6m. Table  1 summarizes the meteorological conditions of the Copenhagen experiment where L is the Obukhov length, h is the height of the convective boundary layer, w * is the convective velocity scale and u * is the friction velocity.
The wind speed profile used in the simulations is described by a power law expressed following the findings of reference Panosfsky & Dutton (1988).
Here u z and u 1 are the horizontal mean wind speed at heights z and z 1 and n is an exponent that is related to the intensity of turbulence (Irwin (1989)). As is possible to see in Irwin (1989)  K z (500,z) K z (1500,z) K z (2500,z) Fig. 1. The vertical eddy diffusivity K z (z, t) for three different travel times (t = 500s, 1500s, 2500s) using run 8 of the Copenhagen experiment. n = 0.1 is valid for a power low wind profile in unstable condition. Moreover, U.S.EPA suggests for rural terrain (as default values used in regulatory models) to use n = 0.15 for neutral condition (class D) and n = 0.1 for stability class C (moderately unstable condition). In Figure 1 we present a plot of K z (z, t) for three different travel times (t = 500s, 1500s, 2500s) using run 8 of the Copenhagen experiment.
In order to exclude differences due to numerical uncertainties we define the numerical accuracy 10 −4 of our simulations determining the suitable number of terms of the solution series. As an eye-guide we report in table 2 on the numerical convergence of the results, considering successively one, two, three and four terms in the solution series. One observes that the desired accuracy, for the solved problem solved is attained including only four terms in the truncated series, which is valid for all distances considered. Once the number of terms in the series solution is determined numerical comparisons of the 3D-GILTT results against experimental data may be performed and are presented in table 3. Figure 2 shows the scatter plot of the centreline ground-level observed concentrations versus the simulated the 3D-GILTT model predictions, normalized by the emission rate and using two points in the time Gaussian Quadrature inversion (Moreira et al. (2006a), Stroud & Secrest (1966)). In the scatter diagram analysis, the closer the data are from the bisector line, the better are the results. The lateral lines indicate a factor of two (FA2), meaning if all the obtained data are between these lines FA2 equals to 1 (the maximum value for stochastically distributed data). From the scatter diagram (figure 2 one observes the fairly good simulation of dispersion data by the 3D-GILTT model, even for the two-point Gaussian quadrature scheme considered here. To perform statistical comparisons between GILTT results against Copenhagen experimental data we consider the set of statistical indices described by Hanna (1989) and defined in by • FA2 = fraction of data (%, normalized to 1) for 0, 5 ≤   Table 4. Statistical comparison between 3D-GILTT model results and the Copenhagen data set, changing the number of terms in equation (18).
where the subscripts o and p refer to observed and predicted quantities, respectively, and the bar indicates an averaged value. The best results are expected to have values near zero for the indices NMSE, FB and FS, and near 1 in the indices COR and FA2. Table 4 shows the findings of the statistical indices that show a fairly good agreement between the 3D-GILTT predictions and the experimental data. Moreover, the splitting proposed for the eddy diffusivity coefficient as a sum of the averaged eddy diffusivity coefficient plus time variation, appears to be a valid assumption, since we got compact convergence of the solution, in the sense that we attained results with accuracy of 10 −4 with only a few terms in the solution series for all the distances considered.

Conclusion
In the present contribution we focused on an analytical description of pollution related phenomena in a micro-scale that allows to simulate dispersion in an computationally efficient procedure. The reason why adopting an analytical procedure instead of using the nowadays available computing power resides in the fact that once an analytical solution to a mathematical model is found one can claim that the problem has been solved. We provided a closed form solution that may be tailored for numerical applications such as to reproduce the solution within a prescribed precision. As a consequence the error analysis reduces to model validation only, in comparison to numerical approaches where in general it is not straight forward to disentangle model errors from numerical ones.
Our starting point, i.e. the mathematical model, is the advection-diffusion equation, which we solved in 3 ⊕ 1 space-time dimensions for a general eddy diffusivity. The closed form solution is obtained using a hybrid approach by spectral theory together with integral transforms (in the present case the Laplace transform) and the decomposition method. The general formalism was simplified in order to attend the meteorological situation of the Copenhagen experiment. By comparison the present approach was found to yield an acceptable solution for the time dependent three dimensional advection-diffusion equation and moreover predicted tracer concentrations closer to observed values compared to other approaches from the literature. Although K-closure is known to have its limitations comparison of measurements and theoretical predictions corresponded on a satisfactory level and thus supported the usage of such an approach for micro-scale dispersion phenomena. Note, that the dispersion of the experimental data in comparison to the predictions might suggest a considerable discrepancy between theory and experiment, but it is worth mentioning that the measurements are a unique sample of a distribution around an average value, whereas the prediction of an average value is evaluated from a deterministic equation, where the stochastic character is hidden in the turbulence closure hypothesis, so that a spread of data along the bisector is to be expected.
A data dispersion due to numerical uncertainties may be excluded using convergence criteria to control the numerical precision. The quality of the solution is controlled by a genuine mathematical convergence criterion. Note, that for the t and x coordinate the Laplace inversion considers only bi-Lipschitz functions, which defines then a unique relation between the original function and its Laplace-transform. This makes the transform procedure manifest exact and the only numerical error comes from truncation, which is determined from the Sturm-Liouville problem. In order to determine the truncation index of the solution series we introduced a carbon-copy of the Cardinal theorem of interpolation theory.
Recalling, that the structure of the pollutant concentration is essentially determined by the mean wind velocityŪ and the eddy diffusivity K, means that the quotient of norms ̟ = ||K|| ||Ū|| defines a length scale for which the pollutant concentration is almost homogeneous. Thus one may conclude that with decreasing length ( ̟ m and m an increasing integer number) variations in the solution become spurious. Upon interpreting ̟ −1 as a sampling density, one may now employ the Cardinal Theorem of Interpolation Theory (Torres (1991)) in order to find the truncation that leaves the analytical solution almost exact, i.e. introduces only functions that vary significantly in length scales beyond the mentioned limit.
The square integrable function χ = rc d td xd η ∈ L 2 (η = y or z) with spectrum {λ i } which is bounded by m̟ −1 has an exact solution for a finite expansion. This statement expresses the Cardinal Theorem of Interpolation Theory for our problem. Since the cut-off defines some sort of sampling density, its introduction is an approximation and is related to convergence of the approach and Parseval's theorem may be used to estimate the error. In order to keep the solution error within a prescribed error, the expansion in the region of interest has to contain n + 1 terms, with n = int mL y,z 2π̟ + 1 2 . For the bounded spectrum and according to the theorem the solution is then exact. In our approximation, if m is properly chosen such that the cut-off part of the spectrum is negligible, then the found solution is almost exact.
Further, the Cauchy-Kowalewski theorem (Courant & Hilbert (1989)) guarantees that the proposed solution is a valid solution of the discussed problem, since this problem is a special case of the afore mentioned theorem, so that existence and uniqueness are guaranteed. It remains to justify convergence of the decomposition method. In general convergence by the decomposition method is not guaranteed, so that the solution shall be tested by an appropriate criterion. Since standard convergence criteria do not apply in a straight forward manner for the present case, we resort to a method which is based on the reasoning of Lyapunov (Boichenko et al. (2005)). While Lyapunov introduced this conception in order to test the influence of variations of the initial condition on the solution, we use a similar procedure to test the stability of convergence while starting from an approximate (initial) solution R 0 (the seed of the recursive scheme). Let |δZ n | = ∑ ∞ i=n+1 R i be the maximum deviation of the correct from the approximate solution Γ n = ∑ n i=0 R i , where · signifies the maximum norm. Then strong convergence occurs if there exists an n 0 such that the sign of λ is negative for all n ≥ n 0 . Here, λ = 1 Γ n log |δZ n | |δZ 0 | . Concluding, analytical solutions of equations are of fundamental importance in understanding and describing physical phenomena, since they might take into account all the parameters of a problem, and investigate their influence. Moreover, when using models, while they are rather sophisticated instruments that ultimately reflect the current state of knowledge on turbulent transport in the atmosphere, the results they provide are subject to a considerable margin of error. This is due to various factors, including in particular the uncertainty of the intrinsic variability of the atmosphere. Models, in fact, provide values expressed as an average, i.e., a mean value obtained by the repeated performance of many experiments, while the measured concentrations are a single value of the sample to which the ensemble average provided by models refer. This is a general characteristic of the theory of atmospheric turbulence and is a consequence of the statistical approach used in attempting to parametrise the chaotic character of the measured data. An analytical solution can be useful in evaluating the performances of numerical models (that solve numerically the advection-diffusion equation) that could compare their results, not only against experimental data but, in an easier way, with the solution itself in order to check numerical errors without the uncertainties presented above. Finally, the program of providing analytical solutions for realistic physical problems, leads us to future problems with different closure hypothesis considering full space-time dependence in the resulting dynamical equation, which we will also approach by the proposed methodology.