Different Numerical Inversion Algorithms of the Laplace Transform for the Solution of the Advection-Diffusion Equation with Non-local Closure in Air Pollution Modeling ^†

ABSTRACT

In this paper, a three-dimensional solution of the steady-state dvection-diffusion equation is obtained applying the so-called generalized integral advection-diffusion multilayer technique, considered non-local closure for turbulent flow. Two different parametrizations were considered for the countergradient coefficient and three different methods of numerical inversion for inverse Laplace transform. The results were compared with the experimental data of Copenhagen experiment by an evaluation of statistical indexes to analyze the solution of the equation through the methods of numerical inversion. Different parametrizations for the vertical turbulent eddy diffusivity and wind profile were utilized. The results show a good agreement with the experiment and the methods of numerical inversion for inverse Laplace transform show almost the same accuracy.

Keywords:
Non-local closure; numerical inversion; pollutant dispersion

RESUMO

Neste trabalho apresenta-se a resolução da equação de advecção-difusão tridimensional estacionária obtida através da técnica GIADMT (Generalized Integral Advection Diffusion Multilayer Technique), considerando o fechamento não-local linear para o fluxo turbulento. Foram consideradas duas parametrizações diferentes para o coeficiente do termo do contragradiente e utilizados três diferentes métodos de inversão numérica para a transformada inversa de Laplace. Comparou-se os resultados com os dados medidos no experimento de Copenhagen através de uma avaliação dos índices estatísticos a fim de comparar a solução da equação através dos métodos de inversão numérica. Ainda, foram utilizados diferentes parametrizações para o coeficiente de difusão turbulento vertical e o perfil do vento. Os resultados apresentaram uma boa concordância com o experimento e os métodos de inversão numérica para a transformada de Laplace apresentaram preaticamente a mesma precisão, sendo que o método baseado na série de Fourier foi o mais acurado dos três algoritmos.Por outro lado, o método de fixed-Talbot foi o que mostrou o melhor desempenho do ponto de vista computacional.

Palavras-chave:
Fechamento não-local; inversão numérica; dispersão de poluentes

1 INTRODUCTION

In the framework of pollutant dispersion modeling in the low atmosphere, Reynolds decomposition applied to the equation resulting from the law of conservation of mass leads to an equation exhibiting the so-called closure problem. In particular, it is necessary to provide constitutive relations linking pollutant concentration to the turbulent fluxes in the atmosphere. The traditional approach is based on the K-theory, or gradient-transport theory, which states, in analogy to molecular diffusion, that the fluxes are proportional to the gradient of the mean pollutant concentration, providing the so-called local Fickian closure³². However, such an approach does not take into account the non-homogeneous character of turbulence in the planetary boundary layer (PBL), which occurs when convective movements are dominant. In fact, K-theory is known to fail in the presence of large-size eddies in the upper portion of the PBL, even though it has been commonly used in models that attempt to describe neutral-to-unstable atmospheric conditions³². Moreover, the derivation of the K-theory based on Prandtl’s mixing-length theory²⁶ is valid only for statically neutral situations³². Various approaches were developed to overcome such a situation. For instance, the nonlocal theories of spectral diffusivity³ and transilient turbulence³¹, which, although supported by some experimental evidence, produce results that differ from those of the statistical dispersion theory³². An alternative consist of considering the so-called nonlocal closure representing the counter-gradient fluxes to indicate the presence of scale eddies or nonlocal fluxes in the PBL^{10), (11), (12}. In this context, one approach considers a Taylor expansion of the turbulent fluxes in the vertical direction^{33), (35}. Here, we consider a linear closure for the counter-gradient term which is consistent with various experiment-based parametrizations of the counter-gradient fluxes^{9), (10), (11), (12), (27}.

From the methodological point of view, the boundary/initial-value problems for advection-diffusion equations with variable coefficients that model pollutant dispersion in the atmosphere are usually solved via integral transform-based methods^{5), (6), (19), (20), (21), (34}. Among those methods, the advection-diffusion multilayer method (ADMM)^{19), (20} provides accurate analytical solutions. In fact, it was reported that the ADMM produces estimations of pollutant concentrations which are as accurate as other integral transform-based methods such as the generalized integral Laplace transform technique (GILTT)^{21), (34} but with remarkably less computational cost²². This last feature is essential in real-life situations such as industrial/natural disasters which require swift (ideally real-time) and accurate estimations of the ground-level distribution and concentration of the pollutants escaped to the atmosphere.

The ADMM is based on the stepwise approximation of the continuous wind velocity profile and eddy diffusivity coefficients in the vertical direction. Such an approximation is the local average of the variable coefficients over each sublayer. Then, the original problem with continuous coefficients is approximated by a problem with piecewise-constant coefficients and interlayer continuity conditions. The approximate problem is then analytically solved in the Laplace space, and the approximation of the solution of the original problem is obtained by applying the inverse Laplace transform to the analytical solution of the approximate problem in the Laplace space. For further details, see²⁰.

On the other hand, the GILTT relies on the Fourier method to produce an truncated eigenfunction series expansion of the solution. As usual, the eigenfunctions are obtained from the related Sturm-Liouville problem resulting via separation of variables. The coefficients of the series solution are obtained by solving an auxiliary problem for a matrix ordinary differential equation via the Laplace transform. This problem is obtained by using the orthogonality property of the eigenfunctions. Finally, the inverse Laplace transform is applied to this approximate solution in the Laplace space in order to obtain the approximation of the solution of the original problem. For further details, see²¹.

Here, we use the generalized integral advection-diffusion multilayer technique (GIADMT) to solve a steady-state three-dimensional advection-diffusion model with nonlocal counter-gradient closure. In particular, the GIADMT applies the GILTT in the crosswind lateral direction to the approximate problem resulting from the application of the ADMM. The resulting multilayer problem for the coefficients of the eigenfunction series expansion of the approximate solution is then solved in the Laplace space. Again, the approximation of the solution of the original problem follows by applying the inverse Laplace transform. For further details, see^{6), (7}.

As the solutions in the Laplace space provided by the three aforementioned methods are very complex, it is usually necessary to take a computational approach to the inversion of the Laplace transform. In models adopting the Fickian closure, we observed²⁸ that, when the solution is sought via the ADMM, the fixed-Talbot inversion algorithm¹ provides more accurate results than the Gaussian quadrature approach³⁰ and than a Fourier series-based method⁸. However, to the best of our knowledge, in models adopting non-Fickian closures, it has not been established which inversion algorithm of the Laplace transform is the most precise.

The objective of this work is twofold: to solve a steady-state three-dimensional advection-diffusion model with non-Fickian counter-gradient closure via the GIADMT, and to establish which inversion algorithm of the Laplace transform is the most accurate in this context. Here, we consider three inversion algorithms, namely, the Gaussian quadrature method³⁰, the fixed-Talbot method¹, and a Fourier series-based method⁸. Also, different parametrizations for the eddy diffusivities, the mean wind velocity profile, and the counter-gradient flux are considered in the numerical simulations, and the corresponding results are compared to the Copenhagen experiment dataset¹⁵ via statistical analysis¹⁷ in order to establish which parametrizations and inversion algorithm provide the most accurate results.

This work is organized as follows: section 2 is devoted to the derivation of the model; in section 3, the formalism of the GIADMT applied to the model is presented including the three inversion algorithms of the Laplace transform; the statistical validation of the model is described for various parametrizations of the coefficients and the the three inversion algorithms in section 4; the accuracy of the computational results is discussed in section 5; and concluding remarks are presented in section 6.

2 PROBLEM STATEMENT

Following^{32), (36}, the Eulerian approach (i.e. fixed reference system) to air pollution modeling is based on the law of conservation of mass for one pollutant species with concentration c(x,y,z,t):

$\frac{\partial c}{\partial t}+u\cdot \nabla c-D\Delta c=S,$ (2.1)

where (x,y,z) ∊ ℝ₊× [-L _y , L _y ] ×[0, h] are the spatial coordinates with L _y , h > 0, t is the temporal coordinate, u is the wind velocity vector field, D is the molecular diffusivity, and S is the pollutant source.

In order to simplify equation (2.1), Reynolds decompositions are assumed for both the wind velocity field and the pollutant concentration, so u = ū +u' and $c-\overline{c}+c\text{'}$, where $\overline{(\cdot )}$ and (·)' represent the mean and turbulent (i.e. fluctuating) parts, respectively. Such decompositions are justified by the existence of the so-called spectral gap, which is the lack of variation at temporal or spatial mesoscales and separates macroscalar mean motions from microscalar turbulent ones (for further details, see sections 2.2 and 2.3 of³²). In addition, it is often assumed^{32), (36} that an ergodic hypothesis² is satisfied by the turbulence (i.e. it is homogeneous and stationary, both statistically), so $〈\overline{(\cdot )}〉=\overline{(\cdot )}$ and 〈(·)'〉=0, where 〈·〉 denotes the Reynolds average operator³². Also, turbulent motions smaller that the mesoscale, as the ones considered here, generally satisfy the conditions⁴ for the so-called incompressibility approximation^{16), (32}, which produces the so-called continuity equation for turbulent fluctuations ∇·u' = 0, implying that u'·∇c' = ∇·(c'u'). With such considerations, and by applying the average operator to equation (2.1), it follows that

$\frac{\partial \overline{c}}{\partial t}+\overline{u}\cdot \nabla \overline{c}+\nabla \cdot 〈c\text{'}u\text{'}〉-D\Delta \overline{c}=〈S〉,$ (2.2)

where 〈c'u'〉 represents the turbulent atmospheric diffusion eddies. Also, as the dispersion effects of molecular diffusion are several orders of magnitude smaller than the ones corresponding to the turbulent diffusion eddies, it is possible to neglect term $D\Delta \overline{c}$ in equation (2.2)^32),(36.

On the other hand, observe that term 〈c'u'〉 introduces three new unknowns, so equation (2.2) needs closure. The traditional approach is the so-called K-theory, or gradient-transport theory, which relies on the so-called Fickian (or first-order local) closure, i.e., 〈c'u'〉 = $-K\nabla \overline{c}$, where K is the second-rank tensor field of turbulent diffusion^{32), (36}. However, when considering point sources in unstable atmospheric conditions, the Fickian closure presents major limitations, as its mixing-length derivation is valid only for statically neutral situations³². In order to overcome such a situation, we consider the nonlocal counter-gradient closure 〈c'u'〉 = $-K\left(\nabla \overline{c}-\gamma \right)$, where γ is the so-called counter-gradient vector field^{10), (11), (12}. With such considerations, and by assuming that the pollutant is nonreactive, so 〈S〉 = S, equation (2.2) becomes

$\frac{\partial \overline{c}}{\partial t}+\overline{u}\cdot \nabla \overline{c}-\nabla \cdot \left(K\left(\nabla \overline{c}-\gamma \right)\right)=S.$ (2.3)

Further simplifications of equation (2.3) can be considered. For instance, K is assumed to be a diagonal tensor with nonzero components K _x , K _y and K _z³⁶. Also, as the dominant convective movements occur in the upward direction, the counter-gradient field γ is usually taken to be aligned to the z-axis, so γ=(0, 0, γ)¹¹. In addition, it is considered that the x-axis is aligned with the wind direction, so that ū =(ū,0,0), and, in consequence, the turbulent diffusion along the x-axis is negligible in comparison to the corresponding advective transport:

$\left|\overline{u}\frac{\partial \overline{c}}{\partial x}\right|\gg \left|\frac{\partial }{\partial x}\left(K_x\frac{\partial \overline{c}}{\partial x}\right)\right|,$ (2.4)

where K _x is the turbulent diffusivity along the x-axis. Also, in presence of a single pollutant source in steady-state emission regime and atmospheric conditions, we have that $\partial \overline{c}/\partial t=0$ and the source term S is treated as a boundary condition. With such considerations, equation (2.3) simplifies to

$\overline{u}\frac{\partial \overline{c}}{\partial x}-\frac{\partial }{\partial y}\left(K_y\frac{\partial \overline{c}}{\partial y}\right)-\frac{\partial }{\partial z}\left(K_z\left(\frac{\partial \overline{c}}{\partial z}-\gamma \right)\right)=\mathrm{0,}$ (2.5)

where K _y and K _z are the turbulent diffusivities along the y- and z-axes, respectively. Equation (2.5) is completed with the boundary condition accounting for the pollutant source

$\overline{u}\overline{c}{|}_{x=0}=Q\delta \left(z-H_s\right)\delta \left(y\right),$ (2.6)

where Q is the emission rate of the pollutant source, which is located at point (0, 0, H _s ), and δ (·) is Dirac’s delta function, and total reflexion conditions³⁶:

$K_y\frac{\partial \overline{c}}{\partial y}{|}_{y=\pm L_y}=\mathrm{0,}$ (2.7)

and

$K_z\left(\frac{\partial \overline{c}}{\partial z}-\gamma \right){|}_{z=\mathrm{0,}h}=0.$ (2. 8)

Finally, observe that the counter-gradient term γ is also unknown, so we propose the linear closure $\gamma =\beta \overline{c}$, which is consistent with various experiment-based parametrizations of γ^{9),(10), (11), (12), (27}. Then, equations (2.5) and (2.8) become, respectively,

$\overline{u}\frac{\partial \overline{c}}{\partial x}-\frac{\partial }{\partial y}\left(K_y\frac{\partial \overline{c}}{\partial y}\right)-\frac{\partial }{\partial z}\left(K_z\frac{\partial \overline{c}}{\partial z}\right)+\frac{\partial }{\partial z}\left(K_z\beta \overline{c}\right)=\mathrm{0,}$ (2.9)

and

$K_z\left(\frac{\partial \overline{c}}{\partial z}-\beta \overline{c}\right){|}_{z=\mathrm{0,}h}=0.$ (2.10)

3 MODEL SOLUTION VIA GIADMT

The solution of problem (2.5)-(2.10) is sought via the GIADMT, which applies the GILTT in variable y to the multilayer problem resulting from the application of the ADMM. This problem for the coefficients of the eigenfunction series expansion of the approximate solution is then solved in the Laplace space. Finally, the approximation of the solution of the original problem follows by applying an inversion algorithm of the Laplace transform.

Following the formalism of the GILTT, the mean pollutant concentration $\overline{c}\left(x,y,z\right)$ is sought as a Fourier series in terms of the eigenfunctions ψ _j (y), where j is the order of the corresponding eigenvalue λ _j , that is,

$\overline{c}\left(x,y,z\right)={\displaystyle \sum _{j=0}^{\infty }\frac{C_j\left(x,z\right){\psi }_j\left(y\right)}{\Vert {\psi }_j\Vert }},\begin{array}{c}\end{array}{\Vert {\psi }_j\Vert }^2={\displaystyle {\int }_{-L_y}^{L_y}{\psi }_j^2\left(y\right)dy.}$ (3.1)

The eigenvalues λ _j and the corresponding eigenfunctions ψ_j (y) in (3.1) are obtained by solving the Sturm-Liouville problem

$\frac{d^2{\psi }_j}{d{y}^2}+{\lambda }_j^2{\psi }_j\left(y\right)=\mathrm{0,}\begin{array}{c}\end{array}\frac{d{\psi }_j}{dy}{|}_{y=\pm L_y}=\mathrm{0,}$ (3.2)

which yields ψ_j (y)= cos(λ _j y), where λ _j = jπ/L _y , so $\Vert {\psi }_j\left(y\right)\Vert =\sqrt{L_y}$.

The coefficients C _j (x, z) of the Fourier series expansion of the mean pollutant concentration $\overline{c}\left(x,y,z\right)$ are obtained by solving the problems resulting from the substitution of (3.1) into (2.9), (2.10) and (2.6), and considering the orthogonality of the eigenfunctions:

$\overline{u}\frac{\partial C_j}{\partial x}+K_y{\lambda }_j^2{C}_j\left(x,z\right)-\frac{\partial }{\partial z}\left(K_z\frac{\partial C_j}{\partial z}\right)+\frac{\partial }{\partial z}\left(K_z\beta C_j\left(x,z\right)\right)=\mathrm{0,}$ (3.3)

$K_z\left(\frac{\partial C_j}{\partial z}-\beta C_j\left(x,z\right)\right){|}_{z=\mathrm{0,}h}=\mathrm{0,}$ (3.4)

$\overline{u}{C}_j\left(\mathrm{0,}z\right)=\frac{Q}{\sqrt{L_y}}\delta \left(z-H_s\right).$ (3.5)

Now, we follow the formalism of the ADMM to find the solution of problem (3.3)-(3.5), that is, divide the height h of the PBL into N sublayers and then take the local average of coefficients K _y , K _z , ū and β in direction z over each sublayer. Explicitly, let ${\left\{z_n\right\}}_{n=\overline{\mathrm{0,}N}}\subset \left[\mathrm{0,}h\right]$ be a partition of the PBL. In each sublayer $\left(z_{n-1},z_n\right),n=\overline{\mathrm{1,}N}$, of thickness ∆z _n = z _n -z _n-1 , consider the following stepwise approximations of ū(z), K _τ (z), τ=y, z, and β(z) given by

${\overline{u}}_n=\frac{1}{\Delta z_n}{\displaystyle {\int }_{z_{n-1}}^{z_n}\overline{u}\left(z\right)dz,}{K}_{{\tau }_n}=\frac{1}{\Delta z_n}{\displaystyle {\int }_{z_{n-1}}^{z_n}{K}_{\tau }\left(z\right)dz,{\beta }_n=\frac{1}{\Delta z_n}{\displaystyle {\int }_{z_{n-1}}^{z_n}\beta \left(z\right)dz.}}$

Let C _jn (x, z)= Cj(x, z) for $\left(x,z\right)\in {ℝ}_{+}^{*}\times \left(z_{n-1},z_n\right),n=\overline{\mathrm{1,}N}$. Then, problem (3.3)-(3.5) is approximated by the following ADMM problem:

${\overline{u}}_n\frac{\partial C_{jn}}{\partial x}+K_{y_n}{\lambda }_j^2{C}_{jn}\left(x,z\right)=K_{z_n}\frac{{\partial }^2{C}_{jn}}{\partial z^2}-K_{z_n}{\beta }_n\frac{\partial C_{jn}}{\partial z},z\in \left(z_{n-1},z_n\right),$ (3.6)

$C_{jn}\left(x,z_n\right)=C_{j\left(n+1\right)}\left(x,z_n\right),$ (3.7)

$K_{z_n}{\left[\frac{\partial C_{jn}}{\partial z}-{\beta }_n{C}_{jn}\left(x,z\right)\right]}_{z=z_n}=K_{z_{n+1}}{\left[\frac{\partial C_{j\left(n+1\right)}}{\partial z}-{\beta }_{n+1}{C}_{j\left(n+1\right)}\left(x,z\right)\right]}_{z=z_n},$ (3.8)

$K_{z_1}{\left[\frac{\partial C_{j1}}{\partial z}-{\beta }_1{C}_{j1}\left(x,z\right)\right]}_{z=0}=K_{z_N}{\left[\frac{\partial C_{jN}}{\partial z}-{\beta }_N{C}_{jN}\left(x,z\right)\right]}_{z=h}=\mathrm{0,}$ (3.9)

$C_{jn}\left(\mathrm{0,}z\right)=\frac{Q}{{\overline{u}}_n\sqrt{L_y}}\delta \left(z-H_s\right){\delta }_{n\overline{n}},z\in \left(z_{n-1},z_n\right),$ (3.10)

in which conditions (3.7) and (3.8) follow from the continuity of the pollutant concentration and the turbulent flux in direction z, respectively, at the partition points z =z _n , $n=\overline{\mathrm{1,}N-1}$. Also, ${\delta }_{n\overline{n}}$ in boundary condition (3.10) is Kronecker’s delta (${\delta }_{n\overline{n}}=1\text{if}n=\overline{n},{\delta }_{n\overline{n}}=0\text{if}n\ne \overline{n}$), and $n=\overline{n}$ indicates the sublayer $\left(z_{\overline{n}-1},z_{\overline{n}}\right)$ in which the pollutant source is located, that is, $H_s\in \left(z_{\overline{n}-1},z_{\overline{n}}\right)$.

By applying the Laplace transform ℒ[·] in direction x to the ADMM problem (3.6)-(3.10), it follows, for each s ∊ ℂ, the ADMM problem in the Laplace space:

$\frac{d^2{\zeta }_{j_n}}{d{z}^2}-{\beta }_n\frac{d{\zeta }_{j_n}}{dz}-\frac{{\overline{u}}_{n}s+K_{y_n}{\lambda }_j^2}{K_{z_n}}{\zeta }_{j_n}\left(s,z\right)=-\frac{Q}{K_{z_n}\sqrt{L_y}}\delta \left(z-H_s\right){\delta }_{n\overline{n}},$ (3.11)

${\zeta }_{j_n}\left(s,z_n\right)={\zeta }_{j_{n+1}}\left(s,z_n\right),$ (3.12)

$K_{z_n}{\left[\frac{\partial {\zeta }_{j_n}}{dz}-{\beta }_n{\zeta }_{j_n}\left(s,z\right)\right]}_{z=z_n}=K_{z_{n+1}}{\left[\frac{d{\zeta }_{j_{n+1}}}{dz}-{\beta }_{n+1}{\zeta }_{j_{n+1}}\left(s,z\right)\right]}_{z=z_n},$ (3.13)

$K_{z_1}{\left[\frac{d{\zeta }_{j_1}}{dz}-{\beta }_1{\zeta }_{j_1}\left(s,z\right)\right]}_{z=0}=K_{z_N}{\left[\frac{d{\zeta }_{j_N}}{dz}-{\beta }_N{\zeta }_{j_N}\left(s,z\right)\right]}_{z=h}=\mathrm{0,}$ (3.14)

where ${\zeta }_{j_n}\left(s,z\right)$ = ℒ [C _jn (x, z)]. Then, for z ∊ (z _n-1 , z _n ), $n=\overline{\mathrm{1,}N}$, the solution of the ADMM problem in the Laplace space (3.11)-(3.14) is

$\begin{array}{l}{\zeta }_{jn}\left(s,z\right)=A_{j_n}{e}^{\left(F_n+R_{jn}\right)z}+B_{j_n}{e}^{\left(F_n-R_{j_n}\right)z}+\\ \begin{array}{ccccccc}& & & & & & \end{array}+\frac{Q}{2R_{j_n}{K}_{z_n}\sqrt{L_y}}\left[e^{\left(F_n-R_{j_n}\right)\left(z-H_s\right)}-e^{\left(F_n+R_{j_n}\right)\left(z-H_s\right)}\right]H\left(z-H_s\right),\end{array}$ (3.15)

where H(·) is Heaviside function,

$F_n=\frac{{\beta }_n}{2},\begin{array}{c}\end{array}{R}_{j_n}=\sqrt{\frac{{\beta }_n^2}{4}+\frac{{\overline{u}}_{n}s+K_{y_n}{\lambda }_j^2}{K_{z_n}}},$ (3.16)

and coefficients $A_{j_n}$ and $B_{j_n}$ in (3.15) are obtained by solving the system of linear equations resulting from the substitution of (3.15) into conditions (3.12)-(3.14).

Then, the solution of the ADMM problem (3.6)-(3.10) follows from the application of the inverse Laplace transform ℒ^-1[·] to (3.15), so the approximate solution to the original problem (2.5)-(2.10) is, formally,

$\overline{c}\left(x,y,z\right)={\displaystyle \sum _{j=0}^{\infty }\frac{\cos \left({\lambda }_{j}y\right)}{\sqrt{L_y}}}{\mathcal{L}}^{-1}\left[{\zeta }_j\right]\left(x,z\right).$ (3.17)

However, note that the complexity of solution (3.15) requires the numerical inversion of the Laplace transform, so the final solution is regarded as semi-analytic. In this work, three inversion algorithms are considered, namely, the Gaussian quadrature method³⁰, the fixed-Talbot method¹, and a Fourier series-based method⁸, in order to establish inversion algorithm provides the most accurate results. The applications of the three algorithms to (3.15) are described next:

• Gaussian quadrature method³⁰:

${\overline{c}}_n\left(x,y,z\right)={\displaystyle \sum _{j=0}^{\infty }\frac{\cos \left({\lambda }_{j}y\right)}{\sqrt{L_y}}}\left\{{\displaystyle \sum _{k=1}^{N_p}\frac{p_k{w}_k}{x}{\zeta }_{jn}\left(\frac{p_k}{x},z\right)}\right\},$ (3.18)

where constants w _k and p _k are, respectively, the weights and roots of the quadrature.

• Fixed-Talbot method¹:

$\begin{array}{l}{c}_n\left(x,y,z\right)={\displaystyle \sum _{j=0}^{\infty }\frac{\cos \left({\lambda }_{j}y\right)}{\sqrt{L_y}}\{\frac{r}{M^{*}}[\frac{1}{2}{\zeta }_{jn}\left(r,z\right)e^{rx}}+\\ \begin{array}{ccccccccc}& & & & & & & & \end{array}+{\displaystyle \sum _{k=1}^{M^{*}-1}Re\left[e^{xS\left({\theta }_k\right)}{\zeta }_{jn}\left(S\left({\theta }_k\right),z\right)\left(1+i\overline{w}\left({\theta }_k\right)\right)\right]}]\},\end{array}$ (3.19)

where r is a parameter with value fixed as r = M*/101x, and S(θ_k )=rθ(cot θ +i), ϖ (θ_k ) = θ_k +(θ_k cot θ_k -1)cot θ_k , θ_k = kπ/M*, -π < θ_k <+ π, and i ² = -1.

• Fourier series-based method⁸:

$\begin{array}{l}{\overline{c}}_n\left(x,y,z\right)={\displaystyle \sum _{j=0}^{\infty }\frac{\cos \left({\lambda }_{j}y\right)}{\sqrt{L_y}}\frac{e^{\alpha x}}{T}[{\displaystyle \sum _{k=1}^{N^{*}}\{Re\left[{\zeta }_{jn}\left(\alpha +\frac{ik\pi }{T},z\right)\cos \left(\frac{k\pi x}{T}\right)\right]-}}\\ \begin{array}{cccccccc}& & & & & & & \end{array}-Im\left[{\zeta }_{jn}\left(\alpha +\frac{ik\pi }{T},z\right)\sin \left(\frac{k\pi x}{T}\right)\right]\}+\frac{1}{2}{\zeta }_{jn}\left(\alpha ,z\right)],\end{array}$ (3.20)

where α and T are free parameter taken here as α =0.0001 and T=55000.

4 MODEL VALIDATION

In order to validate the semi-analytical approach described above, parametrizations must be provided for the eddy diffusivities K _y and K _z , the mean wind velocity profile ū, and the counter-gradient coefficient β. The following parametrizations of are valid for the considered atmospheric convective conditions. For K _z:

• Pleim and Chang ²⁵:

$\frac{K_z}{w_{*}h}=\kappa \frac{z}{h}\left(1-\frac{z}{h}\right),$ (4.1)

where κ = 0.4 is von Kármán constant, and w _* is the convective velocity scale.

• Degrazia et al.¹⁴:

$\frac{K_z}{w_{*}h}=0.22{\left(\frac{z}{h}\right)}^{1/3}{\left(1-\frac{z}{h}\right)}^{1/3}\left[1-e^{-4z/h}-0.0003e^{8z/h}\right],$ (4.2)

• Degrazia et al. ¹³:

$\frac{K_z}{w_{*}h}=\frac{0.09c_w^{1/2}{\psi }^{1/3}{\left(z/h\right)}^{4/3}}{{\left(f_m^{*}\right)}_w^{4/3}}{\displaystyle {\int }_0^{\infty }\frac{\sin \left[\frac{7.84c_w^{1/2}{\psi }^{1/3}{\left(f_m^{*}\right)}_w^{2/3}Xn}{{\left(z/h\right)}^{2/3}}\right]}{{\left(1+n\text{'}\right)}^{5/3}}\frac{dn\text{'}}{n\text{'}},}$ (4.3)

where ψ = 1.26exp(-z/0.8h) is the non-dimensional molecular dissipation rate associated to the plume production, c _v,w = 0.4, $n\text{'}=\widehat{n}\left(1.5z/u{\left(f_m^{*}\right)}_w\right)$ with frequency $\widehat{n}$, X is the non-dimensional time as it is the travel time rate x/u and the convective time scale h/w _*, ${\left(f_m^{*}\right)}_w=z/{\left({\lambda }_m\right)}_w$ is the spectral peak non-dimensional frequency, and ${\left({\lambda }_m\right)}_w=1.8h\left[1-\exp \left(-4z/h\right)-0.0003\exp \left(8z/h\right)\right]$ is the wavelength associated to the maximum of the turbulent vertical spectrum.

Also, for K _y , we use the following parametrization:

• Degrazia et al.¹⁴:

$K_y=\frac{\sqrt{\pi }{\sigma }_{\nu }z}{16{\left(f_m\right)}_{\nu }{q}_{\nu }},$ (4.4)

with

${\sigma }_{\nu }^2=\frac{0.98c_{\nu }}{{\left(f_m\right)}_{\nu }^{2/3}}{\left(\frac{\psi \epsilon }{q_{\nu }}\right)}^{2/3}{\left(\frac{z}{h}\right)}^{2/3}{w}_{*}^2,$ (4.5)

${\psi }_{\epsilon }^{1/3}={\left[{\left(1-\frac{z}{L}\right)}^2{\left(-\frac{z}{L}\right)}^{-2/3}+0.75\right]}^{1/2},$ (4.6)

where u _* is the friction velocity, L is the Monin-Obukhov length, σ _v is the standard deviation of the longitudinal turbulent velocity, (f _m )_v = 0.16 is the lateral wave peak, ψ₍ is the non-dimensional molecular dissipation rate, q _v = 4.16z/h is the stability function, and c _v = 0.4.

For the wind velocity profile ū, we use the two parametrization by Panofsky and Dutton²⁴:

• Power-law profile:

$\overline{u}={\overline{u}}_1{\left(\frac{z}{z_1}\right)}^p,$ (4.7)

where ū and ū ₁ are wind velocities at heights z and z ₁, respectively, whereas p is related to the intensity of the turbulence¹⁸.

• Logarithmic profile:

$\overline{u}=\{\begin{array}{c}\left(u_{*}/\kappa \right)\left[\ln \left(z/z_0\right)-{\Psi }_m\left(z/L\right)\right],\begin{array}{c}\end{array}z\le z_b\\ \left(u_{*}/\kappa \right)\left[\ln \left(z_b/z_0\right)-{\Psi }_m\left(z_b/L\right)\right],\begin{array}{c}\end{array}z\le z_b\end{array},$ (4.8)

where z ₀ is the roughness length of the terrain, z _b = min[|L|, 0.1h] and, with A = [1-16(z/L)]^1/4, Ψ_m is the stability function

${\Psi }_m=\ln \left(\frac{1+A^2}{2}\right)+\ln {\left(\frac{1+A}{2}\right)}^2-2\arctan A+\frac{\pi }{2}.$ (4.9)

For the counter-gradient coefficient β, we use the following two parametrization.

• Cuijpers and Holtslag⁹:

$\beta =\frac{b{w}_{*}^2}{{\sigma }_w^{2}h},$ (4.10)

with constant b and the vertical standard deviation of the turbulent velocity σ _w²⁹:

${\sigma }_w^2=1.8{\left(\frac{z}{h}\right)}^{2/3}{\left(1-\frac{z}{h}\right)}^{2/3}{w}_{*}^2.$ (4.11)

• Roberti et al.²⁷:

$\beta =0.085\frac{q_w}{\Psi h}{\left(\frac{h}{z}\right)}^{2/3},$ (4.12)

with non-dimensional dissipation Ψ=0.913, and q _w is the stability function

$q_w=z{\left[0.594h\left(1-e^{-4z/h}-0.0003e^{8z/h}\right)\right]}^{-1}.$ (4.13)

Model validation is carried out by comparing the simulation results to the Copenhagen experiment observations¹⁵. This experiment consisted of releasing sulfur hexafluoride (SF ₆) from a source of height h = 115m and emission rate Q =100 g/s. Nine experiments were carried out under moderately unstable atmospheric conditions, and data were collected at arches located at 2-6 km from the source. The roughness length is z ₀ = 0.6m. The comparison is performed via statistical analysis using the following indexes¹⁷:

• Normalized mean square error:

$NMSE=\overline{{\left(C_o-C_p\right)}^2}/\overline{C_o{C}_p}\begin{array}{c}\end{array}\left(ideal:NMSE=0\right)$

• Correlation coefficient:

$Cor=\overline{\left[\left(C_o-\overline{C_o}\right)\left(C_p-\overline{C_p}\right)\right]}/{\sigma }_o{\sigma }_p\begin{array}{c}\end{array}\left(ideal:Cor=1\right)$

• Factor of two:

$Fa2=C_p/C_o\in \left[\mathrm{0.5,2}\right]\begin{array}{c}\end{array}\left(ideal:Fa2=1\right)$

• Fractional bias:

$Fb=\left(\overline{C_o}-\overline{C_p}\right)/\left(0.5\left(\overline{C_o}+\overline{C_p}\right)\right)\begin{array}{c}\end{array}\left(ideal:Fb=0\right)$

• Fractional standard deviation:

$Fs=\left({\sigma }_o-{\sigma }_p\right)/\left(0.5\left({\sigma }_o+{\sigma }_p\right)\right)\begin{array}{c}\end{array}\left(ideal:Fs=0\right)$

where subscripts o and p denote the experimentally observed and computationally predicted quantities, respectively, σ is the standard deviation, and the overbar represents the mean value.

5 RESULTS

The numerical methods inversion of the Laplace transform, require parameter choices that have direct impact on the performance, e.g. the number of terms in the Fourier series (N*), the number of terms in the Fixed-Talbot expansion (M*), and the number of points used for the Gaussian quadrature (Np).

In order to determine the choice of such parameters was analysed of the computational cost that each method takes to determine the pollutant concentration given by the equation (3.17). For these tests, was utilized the power wind profile, diffusion coefficient K _z of equation (4.3), term of the counter-gradient of (4.10), micrometeorological parameters of the third arc of experiment 3 of Copenhagen¹⁵; y = 0 and z = 0. The value of the concentration observed in this arc is Co = 3.76.

Table 1 shows the time that each inversion method took to calculate the concentration, the value of the concentration obtained in each calculation and the number of eigenvalues required for GITT, and the criterion of the truncation of the sum of GITT is absolute error less than 10^-2.

Table 1:
Computational cost comparison of each inversion method.

Method	Parameter
FT	M* = 100			M* = 200			M* = 500			M* = 1000
	Time	Cp	Na	Time	Cp	Na	Time	Cp	Na	Time	Cp	Na
	70,67	3,27	31	115,30	3,27	26	648,22	3,27	62	1383,73	3,27	67
FS	N* = 100			N* = 200			N* = 500			N* = 1000
	Time	Cp	Na	Time	Cp	Na	Time	Cp	Na	Time	Cp	Na
	72,86	3,60	31	141,88	3,82	32	342,81	3,72	32	741,24	3,73	32
GQ	Np = 4			Np = 8			Np = 10
	Time	Cp	Na	Time	Cp	Na	Time	Cp	Na
	9,76	5,42	38	11,99	3,82	32	13,11	3,52	31
	Np = 12			Np = 16			Np = 20
	Time	Cp	Na	Time	Cp	Na	Time	Cp	Na
	14,41	3,50	31	17,02	3,64	31	20,09	3,75	31

Analyzing the Table (1), it can be seen that the fixed-Talbot method has equal values of Cp independent of the number of expansion terms, in addition, the calculated value for Cp is very close to Co, but it is a method that takes more time to perform the calculations because it needs many eigenvalues. The higher the value of M*, the more eigenvalues are required. Thus, comparing accuracy/cost, M*=100 is sufficient to calculate the concentration values.

The Fourier series-based method presents small changes in the values of Cp by increasing the number of terms in the series, moreover, the values calculated for Cp are also very close to Co. This method takes less time to perform the calculations than the fixed-Talbot method, because it does not need as many eigenvalues as the fixed-Talbot method. Thus, N* is sufficient to calculate the concentration values.

The Gaussian quadrature method is the fastest method to perform the calculations of Cp, but it is not due to the number of eigenvalues but because the number of points used in the present work is a maximum of twenty. The values calculated for Cp vary according to the value of Np chosen, in addition, not always have values close to Co. It can be seen in Stroud-Secrest³⁰ that the magnitude of the real part of the root of the Gaussian quadrature scheme for the inverse Laplace transform increases with N _p (the order of the approximation). Since the solution to the concentration with the Laplace transform has exponential terms, one can readily observe that from the numerical simulation "overflow" appears for the positive exponential argument and "underflow" for the negative argument when N _p assumes very high values. It is important to note that the calculations were performed on an Intel Core i3, 2.53GHz, 3Gb (RAM) computer. Consequently, to avoid overflow and underflow, the values of N _p were restricted to values around twenty. Thus, Np = 8 is sufficient to calculate the concentration values.

Once the parameters of each inversion method have been determined, the table (2) shows the statistical indexes of the approximate solutions (3.18), (3.19) and (3.20) considering all the parametrizations for the eddy diffusivities K _y (equation (4.4)) and K _z (equations (4.1), (4.2) and (4.3)), the wind velocity profile (equations (4.7) and (4.8)), and the counter-gradient coefficient β with β₁ and β₂ denoting the parametrizations in equations (4.10) and (4.12), respectively. The corresponding computational experiments were performed by truncating the concentration series expansion in equations (3.18), (3.19) and (3.20) after 150 terms.

Table 2:
Statistical indexes.

K z ū	Inversion method	NMSE		Cor		Fa2		Fb		Fs
		β1	β2	β1	β2	β1	β2	β1	β2	β1	β2
1	Gauss	0.34	0.42	0.811	0.814	0.739	0.696	0.31	0.40	0.10	0.17
	Talbot	0.36	0.37	0.820	0.823	0.783	0.739	0.35	0.36	0.15	0.17
	Fourier	0.26	0.32	0.834	0.838	0.826	0.826	0.24	0.32	0.05	0.13
2	Gauss	0.30	0.36	0.808	0.811	0.783	0.739	0.24	0.33	0.02	0.10
	Talbot	0.31	0.32	0.815	0.819	0.783	0.783	0.28	0.29	0.08	0.10
	Fourier	0.24	0.28	0.830	0.834	0.826	0.826	0.18	0.26	-0.02	0.06
3	Gauss	0.31	0.39	0.840	0.846	0.783	0.739	0.31	0.41	0.14	0.22
	Talbot	0.36	0.37	0.841	0.844	0.783	0.783	0.37	0.38	0.21	0.23
	Fourier	0.25	0.31	0.853	0.860	0.826	0.739	0.26	0.35	0.10	0.19
4	Gauss	0.26	0.33	0.834	0.841	0.826	0.696	0.24	0.34	0.05	0.15
	Talbot	0.30	0.31	0.837	0.840	0.826	0.783	0.30	0.32	0.13	0.15
	Fourier	0.22	0.26	0.847	0.855	0.826	0.826	0.19	0.28	0.03	0.12
5	Gauss	0.22	0.25	0.853	0.855	0.826	0.826	0.20	0.27	0.06	0.11
	Talbot	0.22	0.22	0.860	0.863	0.826	0.826	0.22	0.23	0.07	0.08
	Fourier	0.17	0.19	0.871	0.873	0.913	0.870	0.13	0.19	-0.02	0.03
6	Gauss	0.19	0.22	0.851	0.854	0.913	0.913	0.13	0.20	-0.02	0.04
	Talbot	0.20	0.19	0.856	0.859	0.913	0.913	0.16	0.17	0.00	0.01
	Fourier	0.17	0.17	0.867	0.869	0.957	0.913	0.07	0.19	-0.09	-0.04

• logarithmic ū, K z of (4.1); ² power-law ū, K z of (4.1); ³ logarithmic ū, K z of (4.2); ⁴ power-law ū, K z of (4.2); ⁵ logarithmic ū, K z of (4.3); ⁶ power-law ū, K z of (4.3).

It is observed in the Table (2) that the three inversion algorithms studied provide results with similar precision for all parametrizations, since the values of all the statistical indices are close to the ideal values. These three methods are used to calculate the inverse of Laplace in this context of estimating the concentration of pollutants in the atmosphere. However, we can say that the Fourier series-based inversion algorithm produces better results even when comparing the computational effort with the other methods.

Figure 1 shows scatter plots of the observed pollutant concentrations (C _o ) from the Copenhagen experiment versus the predicted concentrations (C _p ) for the power-law wind velocity profile (equation (4.7)), the vertical eddy diffusivity from equation (4.3) and the two counter-gradient coefficients in equations (4.10) and (4.12) via the three inversion algorithms. It is observed that the differences between the inversion algorithms of the Laplace transform are minimal. All the computational experiment showed that the scatter plots in Figure 1 are representative of all the other parametrizations.

Figure 1:
Scatter plots for power-law wind and K _Z of ¹³.

6 CONCLUSIONS

We presented the solution of the steady-state three-dimensional advection-diffusion equation with linear nonlocal closure via the GIADMT with two different parametrizations of the counter-gradient coefficient. Three numerical inversion algorithms of the Laplace transform were evaluated: the Gaussian quadrature method, the fixed-Talbot method, and a Fourier series-based method. The use of non-local closure allowed to model satisfactorily the pollutant concentrations of the Copenhagen experiment independently of the choice of parametrizations and inversion algorithms. The statistical indexes showed that the accuracies of the three algorithms were similar, but the Fourier series-based method was the least expensive from the computational point of view.

ACKNOWLEDGMENTS

Financial support by FAPERGS is gratefully acknowledged.

REFERENCES

• ¹
  J. Abate & P. Valkó. Multi-precision Laplace transform inversion. Int. J. Numer. Methods Eng., 60 (2004), 979-993.
• ²
  M. Beran. Statistical Continuum Theories. John Wiley & Sons, New York, (1968).
• ³
  R. Berkowicz & L. P. Prahm. Generalization of K-theory for turbulent diffusion. Part 1: Spectral turbulent diffusivity concept, J. Appl. Meteor., 18(3) (1976), 266-272.
• ⁴
  J. A. Businger. Equations & Concepts, in: F. T. M. Nieuwstadt, H. van Dop (Eds.) Atmospheric Turbulence and Air Pollution Modelling. Reidel, Dordrecht, (1982), 1-36.
• ⁵
  M. Cassol, S. Wortmann & U. Rizza. Analytical modeling of two-dimensional transient atmospheric pollutant dispersion by double GITT and Laplace Transform techniques, Environ. Model. Soft., 24 (2009), 144-151.
• ⁶
  C. P. Costa, M. T. Vilhena, D. M. Moreira & T. Tirabassi. Semi-analytical solution of the steady three-dimensional advection-diffusion equation in the planetary boundary layer. Atmos. Environ., 40 (2006), 5659-5669.
• ⁷
  C. P. Costa, T. Tirabassi, M. T. Vilhena & D. M. Moreira. A general formulation for pollutant dispersion in the atmosphere. J. Eng. Math., 74(1) (2012), 159-173.
• ⁸
  K. S. Crump. Numerical Inversion of Laplace Transforms Using a Fourier series approximation. J. ACM, 23(1) (1976), 89-96.
• ⁹
  J. W. M. Cuijpers & A. A. M. Holtslag. Impact of skewness and nonlocal effects on scalar and buoyancy fluxes in convective boundary layers, J. Atmos. Sci., 55 (1998), 151-162.
• ¹⁰
  J. W. Deardorff. On the direction and divergence of the small-scale turbulent heat flux, J. Meteor., 18 (1961), 540-548.
• ¹¹
  J. W. Deardorff. The counter-gradient heat flux in the lower atmosphere and in the laboratory, J. Atmosph. Sci., 23 (1966), 503-506.
• ¹²
  J. W. Deardorff. Theoretical expression for the countergradient vertical heat flux, J. Geophys. Res., 77 (1972), 5900-5904.
• ¹³
  G. A. Degrazia, D. M. Moreira & M. T. Vilhena. Derivation of an eddy diffusivity depending on source distance for vertically inhomogeneous turbulence in a convective boundary layer, J. Appl. Meteor. , (2001), 1233-1240.
• ¹⁴
  G. A. Degrazia, U. Rizza, C. Mangia & T. Tirabassi. Validation of a new turbulent parameterization for dispersion models in convective conditions, Boundary-Layer Meteor., 85(2) (1997), 243-254.
• ¹⁵
  S. E. Gryning & E. Lyck. The Copenhagen Tracer Experiments: Reporting of Measurements. Riso National Laboratory, (2002).
• ¹⁶
  S. R. Hanna, G. A. Briggs & R. P. Hosker Jr. Handbook on Atmospheric Diffusion. U. S. Department of Energy, (1982).
• ¹⁷
  S. R. Hanna. Confidence limits for air quality models, as estimated by bootstrap and jackknife resampling methods. Atmos. Environ. , 23 (1989), 1385-1395.
• ¹⁸
  J. S. Irwin. A theoretical variation of the wind profile power-law exponent as a function of surface roughness and stability. Atmos. Environ. , 13(1) (1979), 191-194.
• ¹⁹
  D. M. Moreira, U. Rizza, M. T. Vilhena & A. Goulart. Semi-analytical model for pollution dispersion in the planetary boundary layer. Atmos. Environ. , 39(14) (2005), 2689-2697.
• ²⁰
  D. M. Moreira, M. T. Vilhena, T. Tirabassi, C. P. Costa & B. Bodmann. Simulation of pollutant dispersion in the atmosphere by the Laplace transform: The ADMM approach, Water, Air, Soil Pollut., 77(1) (2006), 411-439.
• ²¹
  D. M. Moreira, M. T. Vilhena, D. Buske & T. Tirabassi. The state-of-art of the GILTT method to simulate pollutant dispersion in the atmosphere, Atmos. Res., 92(1) (2009), 1-17.
• ²²
  D. M. Moreira, M. T. Vilhena, T. Tirabassi, D. Buske & C. P. Costa. Comparison between analytical models to simulate pollutant dispersion in the atmosphere. Int. J. Environ. Waste Manag., 6, No. 3-4 (2010), 327-344.
• ²³
  D. M. Moreira, A. C. Moraes, A. G. Goulart & T. T. A. Albuquerque. A contribution to solve the atmospheric diffusion equation with eddy diffusivity depending on source distance, Atmos. Environ. , 83 (2014), 254-259.
• ²⁴
  A. H. Panofsky & A. J. Dutton. Atmospheric Turbulence, John Wiley & Sons, New York, 1984.
• ²⁵
  J. Pleim & J. Chang. A non-local closure model for vertical mixing in the convective boundary layer, Atmos. Environ. , 26A(6) (1992), 965-981.
• ²⁶
  L. Prandtl. Bericht über Untersuchungen zur ausgebildeten Turbulenz, Z. angew. Math. Mech., 5(2) (1925), 136-139.
• ²⁷
  D. R. Roberti, H. F. Campos Velho & G. A. Degrazia. Identifing counter-gradient term in atmospheric convective boundary layer. Inverse Probl. Sci. Eng., 12(3) (2004), 329-339.
• ²⁸
  K. Rui & C. P. Costa. Comparison of different numerical algorithms for the inverse Laplace transform in the advection-diffusion equation, (2017) (submitted).
• ²⁹
  Z. Sorbjan. Structure of the Atmospheric Boundary Layer. Prentice Hall, (1989).
• ³⁰
  A. H. Stroud & D. Secrest. Gaussian Quadrature Formulas. Prentice Hall, Inc., Englewood Cliffs, N. J., (1966).
• ³¹
  R. B. Stull. Transilient turbulence theory. Part 1: The concept of eddy mixing across finite distances, J. Atmos. Sci., 41(23) (1984), 3351-3367.
• ³²
  R. B. Stull. An Introduction to Boundary Layer Meteorology. Kluwer Academic Publishers, Dordrecht, Holanda, (1988).
• ³³
  H. van Dop & G. Verver. Countergradient transport revisited. J. Atmos. Sci. , 58(15) (2001), 2240-2247.
• ³⁴
  S. Wortmann, M. T. Vilhena, D. M. Moreira & D. Buske. A new analytical approach to simulate the pollutant dispersion in the PBL. Atmos. Environ. , 39(12) (2005), 2171-2178.
• ³⁵
  J. Wyngaard & J. Weil. Transport asymmetry in skewed turbulence. Phys. Fluids A, 3(1) (1991), 155-162.
• ³⁶
  P. Zannetti. Air Pollution Modeling: Theories, Computational Methods and Available Software. Springer, New York, (1990).

• †
  Paper presented at the National Congress of Applied and Computational Mathematics 2016.

Publication Dates

• Publication in this collection
  Jan-Apr 2018

History

• Received
  20 Nov 2016
• Accepted
  18 Oct 2017