%% Antes de processar este arquivo LaTeX (LaTeX2e) deve
%% verificar que o arquivo TEMA.cls estah no mesmo
%% diretorio. O arquivo TEMA.cls pode ser obtido do
%% endereco www.sbmac.org.br/tema.

\documentclass{TEMA}

%\usepackage[brazil]{babel}      % para texto em Português
\usepackage[english]{babel}    % para texto em Inglês

\usepackage[latin1]{inputenc}   % para acentuação em Português
%\input{P-margin.inf}

%\usepackage[dvips]{graphics}
\usepackage{subfigure}
\usepackage{graphicx}
\usepackage{epsfig}

%*********************************************************************************
%comandos adicionados
\usepackage{amsmath,mathrsfs,amsfonts,amssymb,amsthm,bm}
%\usepackage{longtable}
%\usepackage{cite}
%\usepackage[usenames,dvipsnames]{color}
\usepackage{threeparttablex}
\usepackage{placeins}
\usepackage{float}

%\usepackage{enumitem}
%\setlist[itemize]{leftmargin=*}

%*********************************************************************************

\newcommand{\B}{{\tt\symbol{92}}}
\newcommand{\til}{{\tt\symbol{126}}}
\newcommand{\chap}{{\tt\symbol{94}}}
\newcommand{\agud}{{\tt\symbol{13}}}
\newcommand{\crav}{{\tt\symbol{18}}}

\begin{document}

%********************************************************
\title
    {Different Numerical Inversion Algorithms of the Laplace Transform for the Solution of the Advection-Diffusion Equation with Non-local Closure in Air Pollution Modeling\thanks{Financial support by FAPERGS is gratefully acknowledged; Paper presented at the National Congress of Applied and Computational Mathematics 2016.}}

\author
    {Autor %C. P. DA COSTA%
    % \thanks{camiladacosta@gmail.com; doctor teacher.}\,,
    % Departamento de Matemática e Estatística, DME - UFPel,
    % Campus Universitário, 354, 96010-900 Pelotas, RS, Brasil
     \\ \\
    % K. RUI%
    % \thanks{karinerui@gmail.com; PhD student.}\,,
    % Programa de Pós-Graduação em Engenharia Mecânica, PROMEC - UFRGS,
    % Rua Sarmento Leite, 425, 90050-170 Porto Alegre, RS, Brasil
		\\ \\
    % L. D PÉREZ-FERNÁNDEZ%
    % \thanks{leslie.fernandez@ufpel.edu.br;  professor Doctor.}\,,
    % Departamento de Matemática e Estatística, DME - UFPel,
    % Campus Universitário, 354, 96010-900 Pelotas, RS, Brasil
    }

\criartitulo

%\runningheads {Costa, Rui e Pérez-Fernández}{Advection-Diffusion Equation with Non-local Closure}
\runningheads {Autor}{Advection-Diffusion Equation with Non-local Closure}

\begin{abstract}
{\bf Abstract}. In this paper, a three-dimensional solution of the steady-state advection-diffusion equation is obtained applying the Generalized Integral Advection Diffusion Multilayer Technique (GIADMT), considering non-local closure for turbulent flow. Two different parameterizations were considering for the countergradient term and 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 indices to analyse the solution of the equation through the methods of numerical inversion. Differents parameterizations 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 same efficacy.

{\bf Keywords}. Non-local closure, numerical inversion, pollutant dispersion.
\end{abstract}

%********************************************************
\newsec{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 closure \cite{stull1988}. 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 order to study such a situation, the so-called non-local closure representing the counter-gradient fluxes is considered to indicate the presence of scale eddies or non-local fluxes in the PBL. 

This approach produces a more complex advection-diffusion equation which, in this work, is solved via the Generalized Integral Advection-Diffusion Multilayer Technique (GIADMT) for the steady-state case\cite{Costaetal2006}. This method relies on the Laplace transform, whose inversion must be carried out numerically in view of the complexity of the solutions in the Laplace space for this type of problems. For this reason, here we consider three algorithms for the numerical inversion of the Laplace transform, namely, the Gaussian quadrature method \cite{stroudesecrest1966}, the fixed-Talbot method \cite{AbateeValko2004}, and a Fourier series-based method \cite{crump1976}. Also, different parametrizations for the eddy diffusivities, the mean wind velocity profile, and the counter-gradient flux are considered in the simulations, and the corresponding results are compared to the Copenhagen experiment dataset \cite{GryningeLyck2002}. 

\newsec{Problem Statement} 

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)$:
\begin{equation}
\frac{\partial c}{\partial t} + \mathbf{u}\cdot\nabla c - D\Delta c = S,
\label{eq:rui1}
\end{equation}
where $(x,y,z)\in\mathbb{R}_+\times[-L_y,L_y]\times[0,h]$ are the spatial coordinates with $L_y,h>0$, $t$ is the temporal coordinate, $\mathbf{u}$ is the wind velocity vector field, $D$ is the molecular diffusivity, and $S$ is the pollutant source. It is assumed that the velocity field can be decomposed as $\mathbf{u} = \bar{\mathbf{u}} + \mathbf{u}'$, where $\bar{\mathbf{u}}$ and $\mathbf{u}'$ are the mean and fluctuating parts, respectively. Also, a similar decomposition is assumed for the concentration: $c = \langle c \rangle + c'$, where the angular brackets denote the so-called ensemble average. Following the ergodic hypothesis, it assumed that $\langle\bar{\mathbf{u}}\rangle = \bar{\mathbf{u}}$ and $\langle\mathbf{u}'\rangle = 0$. With such considerations, equation (\ref{eq:rui1}) becomes
\begin{equation}
\frac{\partial \langle c \rangle}{\partial t} + \bar{\mathbf{u}}\cdot\nabla \langle c \rangle + \nabla\cdot\langle c'\mathbf{u}' \rangle - D\Delta \langle c \rangle = \langle S \rangle,
\label{eq:rui2}
\end{equation}
where $\langle c'\mathbf{u}' \rangle$ represents the turbulent atmospheric diffusion eddies, with dispersion effects several orders of magnitude larger than the corresponding to molecular diffusion, which leads to neglecting the latter. On the other hand, term $\langle c'\mathbf{u}' \rangle$ introduces three new unknowns, so equation (\ref{eq:rui2}) needs closure. The traditional approach is the so-called $K$-theory, or gradient-transport theory, which relies on the so-called Fickian closure, i.e., $\langle c'\mathbf{u}' \rangle = - \mathbf{K}\nabla\langle c \rangle$, where $K$ is the second-rank tensor field of turbulent diffusion. However, when considering point sources in unstable atmospheric conditions, the Fickian closure presents major limitations. So, in order to overcome such a situation, we consider the nonlocal counter-gradient closure given by $\langle c'\mathbf{u}' \rangle = - \mathbf{K}(\nabla\langle c \rangle - \bm{\gamma})$, where $\bm{\gamma}$ is the so-called counter-gradient vector field. With such considerations, and by assuming that the pollutant is nonreactive, so $\langle S \rangle = S$, equation (\ref{eq:rui2}) becomes
\begin{equation}
\frac{\partial \langle c \rangle}{\partial t} + \bar{\mathbf{u}}\cdot\nabla \langle c \rangle - \nabla\cdot(\mathbf{K}(\nabla\langle c \rangle - \bm{\gamma})) = S.
\label{eq:rui3}
\end{equation}
Further simplifications of equation (\ref{eq:rui3}) can be considered. For instance, $\mathbf{K}$ is assumed to be a diagonal tensor with nonzero components $K_x$, $K_y$ and $K_z$. Also, the counter-gradient field $\bm{\gamma}$ is usually taken to be  aligned to the $z$-axis, so $\bm{\gamma} = (0,0,\gamma)$. In addition, it is considered that the $x$-axis is aligned with the wind direction, so that $\bar{\mathbf{u}} = (\bar{u},0,0)$, and, in consequence, the turbulent diffusion along the $x$-axis is negligible in comparison to the corresponding advective transport:
\begin{equation*}
\left|\bar{u}\frac{\partial \langle c \rangle}{\partial x}\right| \gg \left|\frac{\partial}{\partial x}\left(K_x\frac{\partial \langle c \rangle}{\partial x}\right)\right|,
\label{eq:rui4}
\end{equation*}
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, $\partial \langle c \rangle/\partial t = 0$ and the source term $S$ is treated as a boundary condition. With such considerations, equation (\ref{eq:rui3}) simplifies to
\begin{equation}
\bar{u}\frac{\partial \langle c \rangle}{\partial x} - \frac{\partial}{\partial y}\left(K_y\frac{\partial \langle c \rangle}{\partial y}\right) - \frac{\partial}{\partial z}\left(K_z\left(\frac{\partial \langle c \rangle}{\partial z} - \gamma\right)\right) = 0,
\label{eq:rui5}
\end{equation}
where $K_y$ and $K_z$ are the turbulent diffusivities along the $y$- and $z$-axes, respectively. Equation (\ref{eq:rui5}) is completed with the boundary condition accounting for the pollutant source
\begin{equation}
\left.\bar{u}\langle c \rangle\right|_{x = 0} = Q\delta(z - H_s)\delta(y),
\label{eq:rui6}
\end{equation}
where $Q$ is the emission rate of the pollutant source, which is located at point $(0,0,H_s)$, and $\delta(\cdot)$ is Dirac's delta function, and total reflexion conditions given by
\begin{equation}
\left.K_y\frac{\partial \langle c \rangle}{\partial y}\right|_{y = \pm L_y} = 0,
\label{eq:rui7}
\end{equation}
and
\begin{equation}
\left.K_z\left(\frac{\partial \langle c \rangle}{\partial z} - \gamma\right)\right|_{z = 0,h} = 0.
\label{eq:rui8}
\end{equation}
A final simplification is obtained by assuming that $\gamma = \beta\langle c\rangle$, so equations (\ref{eq:rui5}) and (\ref{eq:rui8}) become
\begin{equation}
\bar{u}\frac{\partial \langle c \rangle}{\partial x} - \frac{\partial}{\partial y}\left(K_y\frac{\partial \langle c \rangle}{\partial y}\right) - \frac{\partial}{\partial z}\left(K_z\frac{\partial \langle c \rangle}{\partial z}\right) + \frac{\partial}{\partial z}\left(K_z\beta\langle c\rangle\right) = 0,
\label{eq:rui9}
\end{equation}
and
\begin{equation}
\left.K_z\left(\frac{\partial \langle c \rangle}{\partial z} - \beta\langle c\rangle\right)\right|_{z = 0,h} = 0.
\label{eq:rui10}
\end{equation}
respectively. In what follows, we denote $\langle c\rangle$ by $\overline{c}$.

\newsec{Model Solution via GIADMT}
%\setcounter{equation}{0}


The solution of problem (\ref{eq:rui5})-(\ref{eq:rui10}) is sought via the GIADMT \cite{Costaetal2006}, which consists of solving the problem resulting from the application of the GITT \cite{COTTA1993} via the ADMM \cite{MOREIRAetal2006}.

The GITT is an integral transform technique that combines a series expansion with an integration. An auxiliary Sturm-Liouville problem is solved and the pollutant concentration is expanded in a Fourier series, which in turn is substituted into the advection-diffusion equation. By applying the orthogonality property of the related eigenfunctions, an infinite system of differential equations is obtained and solved numerically. In this work, the two-dimensional equations in the system are solved via the ADMM, which consists of taking piecewise-constant approximations of the variable coefficients. Then, the resulting equations, coupled with continuity conditions, are solved with the aid of the Laplace transform.

Following the formalism of the GITT, the mean pollutant concentration $\overline{c}(x,y,z)$ is sought as a Fourier series in terms of the orthogonal eigenfunctions $\psi_{j}(y)$ in direction $y$, where $j$ is the order of the corresponding eigenvalues $\lambda_j$, that is, 
\begin{equation}\label{eqgitt}
\overline{c}(x,y,z)=\sum_{j=0}^{\infty} \frac{ C_{j}(x,z)\psi_{j}(y)}{\left\|\psi_{j}\right\|},
\end{equation}
where $\left\|\psi_{j}\right\|^2\!=\!\int_{-L_y}^{L_y} \psi_{j}^{2}(y)\mathrm{d}y$.

The eigenvalues $\lambda_j$ and the corresponding eigenfunctions $\psi_{j}(y)$ in (\ref{eqgitt}) are obtained by solving the Sturm-Liouville problem
\begin{equation*}\label{eqprobauxiliar}
\frac{d^2\psi_j}{d y^2}+\lambda_j^2\psi_j(y)=0, \quad \left.\frac{d\psi_j}{d y}\right|_{y=\pm L_y}=0,
\end{equation*}
which yields $\psi_j(y)=\cos(\lambda_jy)$, where $\lambda_j$ are the positive roots of $\sin(\lambda_jL_y)=0$. 

By substituting (\ref{eqgitt}) truncated at $j = J-1$, $J\in\mathbb{N}$, into (\ref{eq:rui9}), (\ref{eq:rui10}) and (\ref{eq:rui6}), and considering the orthogonality of the eigenfunctions, $J$ problems are obtained:

\begin{equation}
\overline u \frac{\partial C_j}{\partial x} + K_y\lambda^2_j C_j(x,z) - \frac{\partial}{\partial z}\left(K_z\frac{\partial C_j}{\partial z}\right) + \frac{\partial}{\partial z}\left(K_z\beta C_j(x,z)\right)=0,\label{eqadgitt}
\end{equation}
\begin{equation}
\left.K_{z}\left(\dfrac{\partial C_{j}}{\partial z}-\beta C_{j}\right) \right|_{z = 0,h} = 0,\label{eqadgittA}
\end{equation}
\begin{equation}
\overline{u}\,C_{j}(0,z) = \frac{Q \psi_{j}(y_0)}{\left\|\psi_{j}\right\|} \delta(z-H_{s}),\label{eqadgittCC} 
\end{equation}          
for $j=0,1,\ldots,J-1$.
 
The solution of problem \eqref{eqadgitt}-(\ref{eqadgittCC}) is approximated via the ADMM (\emph{Advection Diffusion Multilayer Method}), which consists of dividing the height $h$ of the PBL into $N$ sublayers, in which the mean value of the coefficients $K_{\tau}$, $\overline{u}$ and $\beta$ in direction $z$ are considered \cite{Costaetal2006,MOREIRAetal2006}. 

Let $\{z_n\}_{n = \overline{0,N}}\subset[0,h]$ be a partition of the PBL. In each sublayer $(z_{n - 1}, z_n)$, $n = \overline{1,N}$, of thickness $\Delta z_n = z_n - z_{n - 1}$, consider the following stepwise approximations of $u(z)$, $K_\tau(z)$, $\tau = y,z$, and $\beta$ given by
\begin{equation*}
\overline{u}_n = \frac{1}{\Delta z_n}\int_{z_{n-1}}^{z_n}{\overline{u}(z)dz} \; , \;
K_{\tau_n} = \frac{1}{\Delta z_n}\int_{z_{n-1}}^{z_n}{K_\tau(z)dz}\; , \;
\beta_{n} = \frac{1}{\Delta z_n}\int_{z_{n-1}}^{z_{n}} \beta(z)\,dz.
\label{uKconst}
\end{equation*}
Let $C_{j_n}(x,z) = C_j(x,z)$, $(x,z)\in\mathbb{R}^*_+\times(z_{n-1},z_n)$, $n = \overline{1,N}$. Then, the GITT problem (\ref{eqadgitt})-(\ref{eqadgittCC}) is approximated by the following ADMM problem:
\begin{align}
\overline u_n \frac{\partial C_{j_n}}{\partial x} + K_{y_n}\lambda^2_j C_{j_n}(x,z) &= K_{z_n}\frac{\partial^2 C_{j_n}}{\partial z^2} - K_{z_n}\beta_n\frac{\partial C_{j_n}}{\partial z}, \mbox{ } z\in(z_{n-1},z_n)\label{opeqa}\\
C_{j_n}(x,z_n) &= C_{n + 1}(x,z_n)\label{ccont}\\
K_{z_n}\left[\frac{\partial C_{j_n}}{\partial z} - \beta_nC_{j_n}(x,z)\right]_{z = z_n} &= K_{z_{n + 1}}\left[\frac{\partial c_{n + 1}}{\partial z} - \beta_{n+1}C_{j_{n + 1}}(x,z)\right]_{z = z_n}\label{fcont}\\
K_{z_1}\left[\frac{\partial C_{j_1}}{\partial z} - \beta_1C_{j_1}(x,z)\right]_{z = 0}&= K_{z_N}\left[\frac{\partial C_{j_N}}{\partial z} - \beta_NC_{j_N}(x,z)\right]_{z = h} = 0\\
C_{j_n}(0,z) &= \frac{Q \psi_{j}(y_0)}{\overline u_n \left\|\psi_{j}\right\|}\delta(z - H_s)\delta_{n\bar{n}}, \mbox{ } z\in(z_{n - 1},z_{n})\label{ns}
\end{align}
in which conditions (\ref{ccont}) and (\ref{fcont}) 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{1,N - 1}$. Also, $\delta_{n\bar{n}}$ in boundary condition (\ref{ns}) is Kronecker's delta ($\delta_{n\bar{n}} = 1$ if $n = \bar{n}$, $\delta_{n\bar{n}} = 0$ if $n \neq \bar{n}$), and $n = \bar{n}$ indicates the sublayer $(z_{\bar{n} - 1},z_{\bar{n}})$ in which the pollutant source is located, that is, $H_s\in(z_{\bar{n} - 1},z_{\bar{n}})$. 

By applying the Laplace transform $\mathscr{L}[\cdot]$ with respect to variable $x$ to the ADMM problem (\ref{opeqa})-(\ref{ns}), it follows, for each $s\in\mathbb{C}$, the ADMM problem in the Laplace space
\begin{align}
\frac{d^2 \zeta_{j_n}}{d z^2} - \beta_n\frac{d \zeta_{j_n}}{d z} - \frac{u_n s + K_{y_n}\lambda^2_j}{K_{z_n}}&\zeta_{j_n}(s,z)  = -\frac{Q \psi_{j}(y_0)}{K_{z_n} \left\|\psi_{j}\right\|}\delta(z - H_s)\delta_{n\bar{n}} \label{opeql}\\
\zeta_{j_n}(s,z_n) &= \zeta_{n + 1}(s,z_n)\label{ccontl}\\
K_{z_n}\left[\frac{d\zeta_{j_n}}{d z} - \beta_n\zeta_{j_n}(s,z)\right]_{z = z_n} &= K_{z_{n + 1}}\left[\frac{d\zeta_{n + 1}}{d z} - \beta_{n+1}\zeta_{j_{n + 1}}(s,z)\right]_{z = z_n}\label{fcontl}\\
K_{z_1}\left[\frac{d\zeta_{j_1}}{d z} - \beta_1\zeta_{j_1}(s,z)\right]_{z = 0}&= K_{z_N}\left[\frac{d\zeta_{j_N}}{d z} - \beta_N\zeta_{j_N}(s,z)\right]_{z = h} = 0\label{bcl}
\end{align}
where $\zeta_{j_n}(s,z)=\mathscr{L}[C_{j_n}(x,z)]$. Then, for $z\in(z_{n - 1},z_{n})$, $n = \overline{1,N}$, the solution of the ADMM problem in the Laplace space (\ref{opeql})-(\ref{bcl}) is
\begin{multline}
\zeta_{j_n}(s,z) = A_{j_n}e^{(F_{n}+R_{j_n})z}+ B_{j_n}e^{(F_{n}-R_{j_n})z}\\+\frac{Q \psi_{j}(y_0)}{2R_{j_n}K_{z_n}\left\|\psi_{j}\right\|}\left[e^{(F_{n}-R_{j_n})(z-H_s)}+e^{(F_{n}+R_{j_n})(z-H_s)}\right]H(z-H_s) \, ,
\label{lsol}
\end{multline}
for $n = 1,\ldots,N$, where $H(\cdot)$ is Heaviside function,
\begin{equation*}
F_{n} = -\frac{\beta_n}{2}, \quad R_{j_n} = \sqrt{\frac{\beta_n^{2}}{4}+\frac{\overline{u}_{n}s+K_{y_n}\lambda_{j}^{2}}{K_{z_n}}},
\end{equation*}
and coefficients $A_{j_n}$ and $B_{j_n}$ in (\ref{lsol}) are obtained by solving the system of linear equations resulting from the substitution of (\ref{lsol}) into conditions (\ref{ccontl})-(\ref{bcl}). 

Then, the solution of the ADMM problem (\ref{opeqa})-(\ref{ns}) follows from the application of the inverse Laplace transform $\mathscr{L}^{-1}[\cdot]$ to (\ref{lsol}), so the approximate solution to the original problem (\ref{eq:rui5})-(\ref{eq:rui10}) is, formally, 
\begin{equation*}
\overline{c}(x,y,z)\approx\sum_{j=0}^{J - 1} \frac{ \cos(\lambda_{j}y)}{\left\|\psi_{j}\right\|}\mathscr{L}^{-1}[\zeta_j](x,z).
\end{equation*}
However, note that the complexity of solution (\ref{lsol}) 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 \cite{stroudesecrest1966}, the fixed-Talbot method \cite{AbateeValko2004}, and a Fourier series-based method \cite{crump1976}, whose applications to (\ref{lsol}) are described next.

\begin{itemize}
\item Gaussian quadrature method \cite{stroudesecrest1966}:
\begin{equation}\label{eq:3DQG}
\overline{c}_{j_n}(x,y,z) \approx \sum_{j=0}^{J - 1} \frac{\cos(\lambda_jy)}{\left\|\psi_{j}\right\|}\left\{\sum_{k=1}^{N_p} \frac{p_kw_k}{x}\zeta_{j_n}\left(\frac{p_k}{x},z\right) \right\}
\end{equation}
for $n=1,\ldots,N$. Constants $w_k$ and $p_k$ are, respectively, the weights and roots of the quadrature, and $N_p=8$ is the number of points of the quadrature.
\item Fixed-Talbot method \cite{AbateeValko2004}:
\begin{multline}\label{eq:3DFT}
\overline{c}_{j_n}(x,y,z) \approx \sum_{j=0}^{J - 1} \frac{\cos(\lambda_jy)}{\left\|\psi_{j}\right\|}\left\{\frac{r}{M^*}\left[\frac{1}{2} \zeta_{j_n}(r,z)e^{rx}\right.\right.\\
\left.\left.+ \sum_{k=1}^{M^*-1}Re\left[e^{xS(\theta_k)} \zeta_{j_n}(S(\theta_k),z)(1+i\overline{w}(\theta_k))\right]\right]\right\}\, ,
\end{multline}
where $r$ is a parameter with value fixed as $r=2M^{*}/101x$, $M^*=100$, $S(\theta_k)=r\theta(\cot{\theta}+ i)$, $\overline{\omega}(\theta_k)= \theta_k+(\theta_k\cot{\theta_k}-1)\cot{\theta_k}$, $\theta_k=k\pi/M^*$, $-\pi<\theta_k<+\pi$, and $i^2= -1$.
\item Fourier series-based method \cite{crump1976}:
\begin{multline}\label{eq:3DF}
\overline{c}_{j_n}(x,y,z) \approx \sum_{j=0}^{J - 1} \frac{\cos(\lambda_jy)}{\left\|\psi_{j}\right\|} \frac{e^{\gamma x}}{T} \left[\sum_{k=1}^{N^*}\left\{Re\biggl[ \zeta_{j_n}\left(\gamma + \frac{ik\pi}{T},z\right)\cos\left(\frac{k\pi x}{T}\right)\biggr] \right.\right.\\
\left.\left.-Im\left[ \zeta_{j_n}\left(\gamma + \frac{ik\pi}{T},z\right)\sin\left(\frac{k\pi x}{T}\right)\right] \right\} + \frac{1}{2} \zeta_{j_n}(\gamma,z)\right] 
\end{multline}
where $\gamma$ and $T$ are free parameter taken here as $\gamma=0.0001$ and $T=55000$ for $N^{*}=1000$.
\end{itemize}

\newsec{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 $\overline{u}$, and the counter-gradient coefficient $\beta$. The following parametrizations of are valid for the considered atmospheric convective conditions. For $K_z$:
\begin{itemize}
\item Pleim and Chang \cite{PleimeChang1992}: 
\begin{equation*}
\frac{K_{z}}{w_*h}=\kappa \frac{z}{h}\left(1-\frac{z}{h}\right),
\end{equation*}
where $\kappa=0.4$ is von K\'{a}rm\'{a}n constant, and $w_*$ is the convective velocity scale.
\item Degrazia et al. \cite{Degraziaetal1997}: 
\begin{equation*}
\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],
\end{equation*}
\item Degrazia et al. \cite{Degraziaetal2001}: 
\begin{equation*}
\frac{K_{z}}{w_*h}=\frac{0.09c_{w}^{1/2}\psi^{1/3}(z/h)^{4/3}}{(f_{m}^{*})_{w}^{4/3}}\int_{0}^{\infty} \frac{\sin\left[\frac{7.84c_{w}^{1/2}\psi^{1/3}(f_{m}^{*})_{w}^{2/3}Xn}{(z/h)^{2/3}}\right]}{(1+n')^{5/3}} \frac{\mathrm{d}n'}{n'}.
\end{equation*}
where $\psi=1.26\exp(-z/0.8h)$ is the non-dimensional molecular dissipation rate associated to the plume production, $c_{v,w}=0.4$, $n'=\hat{n}(1.5z/u(f_{m}^{*})_{w})$ with frequency $\hat{n}$, $X$ is the non-dimensional time as it is the travel time rate $x/u$ and the convective time scale $h/w_{*}$, $(f_{m}^{*})_{w}=z/(\lambda_{m})_{w}$ is the spectral peak non-dimensional frequency, and $(\lambda_{m})_{w}=1.8h[1-\exp(-4z/h)-0.0003\exp(8z/h)]$ is the wavelength associated to the maximum of the turbulent vertical spectrum.
\end{itemize}

Also, for $K_y$, we use the following parametrization:

\begin{itemize}
\item Degrazia et al. \cite{Degraziaetal1997}: 
\begin{equation*}
K_{y}=\frac{\sqrt{\pi}\sigma_{v}z}{16(f_{m})_{v}q_{v}},
\end{equation*}
with
\begin{align*}
\sigma_{v}^{2}=\frac{0.98c_{v}}{(f_{m})_{v}^{2/3}}\left(\frac{\psi_{\epsilon}}{q_{v}}\right)^{2/3}\left(\frac{z}{h}\right)^{2/3}w_{*}^{2}\, ,
\end{align*}
\begin{align*}
\psi_{\epsilon}^{1/3}=\left[\left(1-\frac{z}{L}\right)^{2}\left(-\frac{z}{L}\right)^{-2/3}+0.75\right]^{1/2}\,,
\end{align*}
where $u_*$ is the friction velocity, $L$ is the Monin-Obukhov length, $\sigma_{v}$ is the standard deviation of the longitudinal turbulent velocity, $(f_{m})_{v}=0.16$ is the lateral wave peak, $\psi_{\epsilon}$ is the non-dimensional molecular dissipation rate, $q_{v}=4.16z/h$ is the stability function, and $c_{v}=0.4$.
\end{itemize}

For the wind velocity profile $\overline{u}$, we use the two parametrization by Panofsky and Dutton \cite{PanofskyeDutton1984}: 
\begin{itemize}
\item Power-law profile: 
\begin{equation*}
\frac{\overline{u_{z}}}{\overline{u_{1}}}=\left(\frac{z}{z_1}\right)^{p},
\end{equation*}
where $\overline{u_z}$ and $\overline{u_1}$ are wind velocities at heights $z$ and $z_1$, respectively, whereas $p$ is related to the intensity of the turbulence \cite{Irwin1979}.
\item Logarithmic profile: 
\begin{equation*}
u=\frac{u_{*}}{\kappa}\left[\ln\left(\frac{z}{z_{0}}\right)-\Psi_{m}\left(\frac{z}{L}\right)\right] \quad \mbox{se} \quad z \leq z_b.
\end{equation*} 
where $z_{0}$ is the roughness length of the terrain, $z_{b}=\min[|L|,0.1h]$ and $\Psi_{m}$ is the stability function
\begin{equation*}
	\Psi_{m} = \ln\left(\frac{1+A^2}{2}\right) + \ln\left(\frac{1+A}{2}\right)^{2} - 2\tan^{-1}{(x)} + \frac{\pi}{2},
\end{equation*} 
with $A = \left[1-15(z/L)\right]^{1/4}$.
\end{itemize}

For the counter-gradient coefficient $\beta$, we use the following two parametrization.
\begin{itemize}
\item Cuijpers and Holtslag \cite{CuijperseHoltslag1998}: 
\begin{equation*}
\beta = \frac{bw_{*}^{2}}{\sigma_{w}^{2}h},
\end{equation*}
where $b$ is a constant and $\sigma_{w}$ is the vertical standard deviation of the turbulent velocity by Sorbjan \cite{Sorbjan1989}:
\begin{equation*}
\sigma_{w}^{2}=1.8\left(\frac{z}{h}\right)^{2/3}\left(1-\frac{z}{h}\right)^{2/3}w_{*}^{2}.
\end{equation*}
\item Roberti \cite{Robertietal2004}: 
\begin{equation*}
\beta = 0.085\frac{q_{w}}{\Psi h}\left(\frac{h}{z}\right)^{2/3},
\end{equation*}
where $\Psi=0.913$ is the non-dimensional dissipation function, and $q_{w}$ is the stability function 
\begin{equation*}
q_{w}=z\left[0.594h\left(1-e^{-4z/h}-0.0003e^{8z/h}\right)\right]^{-1}.
\end{equation*}
\end{itemize}

Model validation is carried out by comparing the simulation results to the Copenhagen experiment observations \cite{GryningeLyck2002}. This experiment consisted of releasing sulfur hexafluoride ($SF_6$) from a source of height $h = 115 \, m$ 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 = 0.6 \, m$. The comparison is performed via statistical analysis using the indexes recomended by Hanna \cite{Hanna1989}, namely: 
\begin{itemize}
\item Normalized mean square error: 

$N\!M\!S\!E =\overline{(C_{o}-C_{p})^{2}}/\overline{C_{o}}\overline{C_{p}} \quad (\mbox{ideal:} \, N\!M\!S\!E=0) $

\item Correlation coefficient: 

$Cor=[\overline{(C_{o}-\overline{C_{o}})(C_{p}-\overline{C_{p}})}]/\sigma_{o}\sigma_{p} \quad (\mbox{ideal:} \, Cor=1)$

\item Factor of two: 

$Fa2=C_{p}/C_{o} \in [0.5,2] \quad (\mbox{ideal:} \, Fa2=1) $

\item Fractional bias: 

$Fb=(\overline{C_{o}}-\overline{C_{p}})/(0.5(\overline{C_{o}}+\overline{C_{p}})) \quad (\mbox{ideal:} \, Fb=0) $

\item Fractional standard deviation: 

$Fs=(\sigma_{o}-\sigma_{p})/(0.5(\sigma_{o}+\sigma_{p})) \quad (\mbox{ideal:} \, Fs=0) $
\end{itemize}
\noindent where subscripts $o$ and $p$ denote the experimentally observed and computationally predicted quantities, respectively, and $\sigma$ is the standard deviation. 

\newsec{Results}

Table \ref{tabela01} shows the performance of the three approximate solutions (\ref{eq:3DQG}), (\ref{eq:3DFT}) and (\ref{eq:3DF}) considering all the parametrizations described above. Such a performance is assessed via the statistical indexes also presented above. Also, $\beta_1$ and $\beta_2$ denote the parametrizations of the counter-gradient coefficient by Cuijpers and Holtslag \cite{CuijperseHoltslag1998} and Roberti et al. \cite{Robertietal2004}, respectively. It is observed that the three inversion algorithms provided results with similar accuracy for all the parametrizations considered here, as the values of the statistical indexes are near their ideal values, whereas the Fourier series-based method is slightly better.

\begin{table}[H]
\begin{threeparttable}[b]
\caption{{\small Statistical indexes}}\label{tabela01}
{\footnotesize
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline $K_z$         & Inversion       & \multicolumn{2}{c|}{$N\!M\!S\!E$} & \multicolumn{2}{c|}{$Cor$} & \multicolumn{2}{c|}{$Fa2$} & \multicolumn{2}{c|}{$Fb$} & \multicolumn{2}{c|}{$Fs$} \\\cline{3-12}
       $\overline u$ & method          & $\beta_{1}$                      & $\beta_{2}$               & $\beta_{1}$               & $\beta_{2}$  & $\beta_{1}$ & $\beta_{2}$                      & $\beta_{1}$               & $\beta_{2}$               & $\beta_{1}$  & $\beta_{2}$ \\
\hline
\hline  & Gauss   & $0.34$ & $0.42$ & $0.811$ & $0.814$ & $0.739$ & $0.696$ & $0.31$  & $0.40$ & $0.10$  & $0.17$ \\ \cline{3-12}
		\tnote{1} & Talbot     & $0.36$ & $0.37$  & $0.820$ & $0.823$ & $0.783$ & $0.739$ & $0.35$ & $0.36$  & $0.15$ & $0.17$ \\ 
		\cline{3-12}
		& Fourier & $0.26$ & $0.32$ & $0.834$ & $0.838$ & $0.826$ & $0.826$ & $0.24$  & $0.32$ & $0.05$  & $0.13$ \\
\hline  & Gauss   & $0.30$ & $0.36$ & $0.808$ & $0.811$ & $0.783$ & $0.739$ & $0.24$  & $0.33$ & $0.02$  & $0.10$ \\ 
        \cline{3-12}
		\tnote{2} & Talbot     & $0.31$ & $0.32$  & $0.815$ & $0.819$ & $0.783$ & $0.783$ & $0.28$ & $0.29$  & $0.08$ & $0.10$ \\ 
		\cline{3-12}
		& Fourier & $0.24$ & $0.28$ & $0.830$ & $0.834$ & $0.826$ & $0.826$ & $0.18$  & $0.26$ & -$0.02$ & $0.06$ \\
\hline  & Gauss   & $0.31$ & $0.39$ & $0.840$ & $0.846$ & $0.783$ & $0.739$ & $0.31$  & $0.41$ & $0.14$  & $0.22$ \\ 
        \cline{3-12}
		\tnote{3} & Talbot     & $0.36$ & $0.37$  & $0.841$ & $0.844$ & $0.783$ & $0.783$ & $0.37$ & $0.38$  & $0.21$ & $0.23$\\  
		\cline{3-12}	
		& Fourier & $0.25$ & $0.31$ & $0.853$ & $0.860$ & $0.826$ & $0.739$ & $0.26$  & $0.35$ & $0.10$  & $0.19$ \\
\hline  & Gauss   & $0.26$ & $0.33$ & $0.834$ & $0.841$ & $0.826$ & $0.696$ & $0.24$  & $0.34$ & $0.05$  & $0.15$ \\  
        \cline{3-12}
		\tnote{4} & Talbot     & $0.30$ & $0.31$  & $0.837$ & $0.840$ & $0.826$ & $0.783$ & $0.30$ & $0.32$  & $0.13$ & $0.15$\\  
		\cline{3-12}
		& Fourier & $0.22$ & $0.26$ & $0.847$ & $0.855$ & $0.826$ & $0.826$ & $0.19$  & $0.28$ & $0.03$  & $0.12$ \\
\hline  & Gauss   & $0.22$ & $0.25$ & $0.853$ & $0.855$ & $0.826$ & $0.826$ & $0.20$  & $0.27$ & $0.06$  & $0.11$ \\   
        \cline{3-12}
		\tnote{5} & Talbot     & $0.22$ & $0.22$  & $0.860$ & $0.863$ & $0.826$ & $0.826$ & $0.22$ & $0.23$  & $0.07$ & $0.08$ \\ 
		\cline{3-12}
		& Fourier & $0.17$ & $0.19$ & $0.871$ & $0.873$ & $0.913$ & $0.870$ & $0.13$  & $0.19$ & -$0.02$ & $0.03$ \\
\hline  & Gauss   & $0.19$ & $0.22$ & $0.851$ & $0.854$ & $0.913$ & $0.913$ & $0.13$  & $0.20$ & -$0.02$ & $0.04$ \\ 
        \cline{3-12}
		\tnote{6} & Talbot     & $0.20$ & $0.19$  & $0.856$ & $0.859$ & $0.913$ & $0.913$ & $0.16$ & $0.17$  & $0.00$ & $0.01$\\ 
		\cline{3-12}
		& Fourier & $0.17$ & $0.17$ & $0.867$ & $0.869$ & $0.957$ & $0.913$ & $0.07$  & $0.19$ & -$0.09$ & -$0.04$ \\
\hline
\end{tabular}
}
\begin{tablenotes}\footnotesize
		\item[1] logarithmic $\overline{u}$, $K_z$ of \cite{PleimeChang1992};
		$^2$ power-law $\overline{u}$, $K_z$ of \cite{PleimeChang1992};
		$^3$ logarithmic $\overline{u}$, $K_z$ of \cite{Degraziaetal1997};\\
		$^4$ power-law $\overline{u}$, $K_z$ of \cite{Degraziaetal1997};
		$^5$ logarithmic $\overline{u}$, $K_z$ of \cite{Degraziaetal2001};
		$^6$ power-law $\overline{u}$, $K_z$ of \cite{Degraziaetal2001}.
\end{tablenotes} 
\end{threeparttable}
\end{table}

Figure \ref{fig1} 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, the vertical eddy diffusivity of Degrazia et. al \cite{Degraziaetal2001} and the two counter-gradient coefficients via the three inversion algorithms. It is observed that the differences between the three choices of numerical inversion of the Laplace transform are minimal.

\begin{figure}[H]
\centering %
\subfigure[Cuijpers and Holtslag (1998)]{\scalebox{1}{\includegraphics[height=5.5cm,width=5.5cm]{espalha_ch.eps}} } %
\hspace{1cm}  %
\subfigure[Roberti et al. (2004)]{\scalebox{1}{\includegraphics[height=5.5cm,width=5.5cm]{espalha_r.eps}} } %
\caption{Scatter plots for power-law wind and $K_z$ of \cite{Degraziaetal2001}.} %
\label{fig1}
\end{figure}
\FloatBarrier

\newsec{Conclusions}

We presented the solution of the steady-state three-dimensional advection-diffusion equation with non-local closure via the GIADMT with two different counter-gradient parametrizations. 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 fixed-Talbot method was the least expensive from the computational point of view.

\begin{abstract}
{\bf 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 (\emph{Generalized Integral Advection Diffusion Multilayer Technique}), considerando o fechamento não-local para o fluxo turbulento. Foram consideradas duas parametrizações diferentes para o termo do contragradiente e utilizados 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. Utilizou-se 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 a mesma eficácia.

{\bf Palavras-chave}. Fechamento não-local, inversão numérica, dispersão de poluentes.
\end{abstract}

\begin{thebibliography}{8}

\bibitem{AbateeValko2004} J. Abate, P. Valk\'{o}, Multi-precision Laplace transform inversion, {\em Int. J. Numer. Methods Eng.}, {\bf 60} (2004), 979--993.

\bibitem{Costaetal2006} 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, {\em Atmos. Environ.}, {\bf 40} (2006), 5659--5669.

\bibitem{COTTA1993} R. M. Cotta, ``Integral transforms in computational heat and fluid flow'', Boca Raton, Florida: CRC Press, 1993.

 
\bibitem{crump1976} K. S. Crump, Numerical Inversion of {L}aplace Transforms Using a {F}ourier series approximation, {\em J. ACM}, {\bf 23}, No. 1 (1976), 89--96.

\bibitem{CuijperseHoltslag1998} J. W. M. Cuijpers, A. A. M. Holtslag, Impact of skewness and nonlocal effects on scalar and buoyancy fluxes in convective boundary layers, {\em J. Atmos. Sci.}, {\bf 55} (1998), 151--162.

\bibitem{Deardorff1972} J. W. Deardorff, Theoretical expression for the countergradient vertical heat flux, {\em J. Geophys. Res.}, {\bf 77} (1972), 5900--5904.

\bibitem{Degraziaetal2001} 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, {\em J. Appl. Meteor.}, (2001), 1233--1240.

\bibitem{Degraziaetal1997} G. A. Degrazia, U. Rizza, C. Mangia, T. Tirabassi, Validation of a new turbulent parameterization for dispersion models in convective conditions, {\em Boundary-Layer Meteor.}, {\bf 85}, No. 2 (1997), 243--254.

\bibitem{GryningeLyck2002} S. E. Gryning, E. Lyck, ``The Copenhagen Tracer Experiments: Reporting of Measurements'', Riso National Laboratory, 2002.

\bibitem{Hanna1989} S. R. Hanna, Confidence limits for air quality models, as estimated by bootstrap and jackknife resampling methods, {\em Atmos. Environ.}, {\bf 23} (1989), 1385--1395.

\bibitem{Irwin1979} J. S. Irwin, A theoretical variation of the wind profile power-law exponent as a function of surface roughness and stability, {\em Atmos. Environ.}, {\bf 13}, No. 1 (1979), 191--194.

\bibitem{MOREIRAetal2006} 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 {A}{D}{M}{M} approach, {\em Water, Air, Soil Pollut.}, {\bf 77}, No. 1-4 (2006), 411--439.

\bibitem{moreira14} 
Moreira, D.M., Moraes, A.C., Goulart, A.G., Albuquerque, T.T.A., 2014. A contribution to solve the atmospheric diffusion equation with eddy diffusivity depending on source distance. Atmospheric Environment 83, 254--259. 

\bibitem{PanofskyeDutton1984} A. H. Panofsky, A. J. Dutton, ``Atmospheric Turbulence'', John Wiley \& Sons, New York, 1984.

\bibitem{PleimeChang1992} J. Pleim, J. Chang, A non-local closure model for vertical mixing in the convective boundary layer, {\em Atmos. Environ.}, {\bf 26A}, No. 6 (1992), 965--981.

\bibitem{Robertietal2004} D. R. Roberti, H. F. Campos Velho, G. A. Degrazia, Identifing counter-gradient term in atmospheric convective boundary layer, {\em Inverse Probl. Sci. Eng.}, {\bf 12}, No. 3 (2004), 329--339.

\bibitem{Sorbjan1989} Z. Sorbjan, ``Structure of the Atmospheric Boundary Layer'', Prentice Hall, 1989.

\bibitem{stroudesecrest1966} A. H. Stroud, D. Secrest,  ``Gaussian Quadrature Formulas'', Prentice Hall, Inc., Englewood Cliffs, N. J., 1966.

\bibitem{stull1988} R. B. Stull,  ``An Introduction to Boundary Layer Meteorology'', Kluwer Academic Publishers, Dordrecht, Holanda, 1988.
\end{thebibliography}


\end{document}
\newpage
$ \  \  $  \thispagestyle{myheadings}  \markboth{      }{   }