Fractional derivative models for atmospheric dispersion of pollutants

TL;DR

This paper explores fractional derivative models to improve the accuracy of atmospheric pollutant dispersion predictions, demonstrating superior performance over traditional and existing variable-coefficient models through real-world comparison.

Contribution

It introduces simple fractional differential equation models for pollutant dispersion and shows they outperform traditional Gaussian and variable-coefficient models in real experiments.

Findings

Fractional models outperform Gaussian models in accuracy.

Fractional models outperform variable-coefficient models.

Better fit to real atmospheric dispersion data.

Abstract

In the present work, we investigate the potential of fractional derivatives to model atmospheric dispersion of pollutants. We propose simple fractional differential equation models for the steady state spatial distribution of concentration of a non-reactive pollutant in Planetary Boundary Layer. We solve these models and we compare the solutions with a real experiment. We found that the fractional derivative models perform far better than the traditional Gaussian model and even better than models found in the literature where it is considered that the diffusion coefficient is a function of the position in order to deal with the anomalous diffusion.

Tables

Table 1. Table 1: Meteorological conditions during the Copenhagen experiment

Exp.	${\overline{u}}_{10}\left(m{s}^{-1}\right)$	$u_{\ast }\left(m{s}^{-1}\right)$	L(m)	${\sigma }_w\left(m{s}^{-1}\right)$	$h$ (m)
1	2.1	0.37	-46	0.83	1980
2	4.9	0.74	-384	1.07	1920
3	2.4	0.39	-108	0.68	1120
4	2.5	0.39	-173	0.47	390
5	3.1	0.46	-577	0.71	820
6	7.2	1.07	-569	1.33	1300
7	4.1	0.65	-136	0.87	1850
8	4.2	0.70	-72	0.72	810
9	5.1	0.77	-382	0.98	2090

Table 2. Table 2: Observed and estimated crosswind-integrated concentrations c y ¯ Q ¯ superscript 𝑐 𝑦 𝑄 \frac{\overline{c^{y}}}{Q} ( 10 − 4 ​ s ​ m − 2 superscript 10 4 𝑠 superscript 𝑚 2 10^{-4}sm^{-2} ) for Copenhagen experiment.

Exp.	Distance (m)	Observed	$\alpha$–Gaussian	Gaussian	$\alpha$–SM	SM
1	1900	6.48	6.32	3.61	7.74	4.62
1	3700	2.31	4.97	2.72	4.20	2.30
2	2100	5.38	4.14	2.47	4.89	3.18
2	4200	2.95	3.27	1.76	2.79	1.58
3	1900	8.20	6.51	4.00	8.91	5.66
3	3700	6.22	5.22	3.73	5.03	2.88
3	5400	4.30	4.66	3.72	3.95	2.21
4	4000	11.7	10.60	10.25	8.21	5.93
5	2100	6.72	5.71	3.98	7.55	4.81
5	4200	5.84	4.70	3.93	4.28	2.43
5	6100	4.97	4.36	3.93	3.27	2.00
6	2000	3.96	2.90	1.72	3.90	2.63
6	4200	2.22	2.27	1.24	2.23	1.28
6	5900	1.83	2.08	1.12	1.79	0.90
7	2000	6.70	4.68	2.77	6.21	4.16
7	4100	3.25	3.65	1.95	3.50	2.03
7	5300	2.23	3.34	1.73	2.79	1.56
8	1900	4.16	5.75	3.51	7.32	4.87
8	3600	2.02	4.72	3.01	4.68	2.74
8	5300	1.52	4.18	2.95	3.39	1.84
9	2100	4.58	3.77	2.26	5.26	3.44
9	4200	3.11	2.99	1.61	3.07	1.74
9	6000	2.59	2.63	1.35	2.24	1.19

Table 3. Table 3: Statistical indices to evaluate the performance of proposed models (Eq. 38 )

Model	Cor	NMSE	FS	FB	FA2
$\alpha$–Gaussian	0.83	0.07	0.30	0.001	0.87
Gaussian	0.82	0.23	0.27	-0.39	0.73
$\alpha$–SM	0.83	0.08	0.17	0.03	0.83
SM	0.80	0.30	0.50	-0.44	0.74

Equations

$${}_a{J}_x^{\alpha }f(x)=\frac{1}{\Gamma (\alpha )}{\int }_a^x\frac{f(u)}{(x-u{)}^{1-\alpha }}du(\alpha \in {\mathbb{R}}_{+})$$

$${}_a{D}_x^{\alpha }f(x)=\frac{1}{\Gamma (n-\alpha )}\frac{d^n}{d{x}^n}{\int }_a^x\frac{f(u)}{(x-u{)}^{1+\alpha -n}}du(n=[\alpha ]+1,\alpha \in {\mathbb{R}}_{+}^{\ast },a\in \mathbb{R}).$$

$${}_a^C{D}_x^{\alpha }f(x)=\frac{1}{\Gamma (n-\alpha )}{\int }_a^x\frac{f^{(n)}(u)}{(x-u{)}^{1+\alpha -n}}du(n=[\alpha ]+1,\alpha \in {\mathbb{R}}_{+}^{\ast },a\in \mathbb{R}),$$

$${}_a{D}_x^{\alpha }(x-a{)}^{\beta }=\frac{\Gamma (\beta +1)}{\Gamma (\beta +1-\alpha )}(x-a{)}^{\beta -\alpha },$$

$${}_a^C{D}_x^{\alpha }(x-a{)}^{\beta }=\frac{\Gamma (\beta +1)}{\Gamma (\beta +1-\alpha )}(x-a{)}^{\beta -\alpha }.$$

$$\frac{{\partial }^{\alpha }}{\partial x_i^{\alpha }}f(x_1,...,x_d)=\frac{1}{\Gamma (n-\alpha )}\frac{{\partial }^n}{\partial x_i^n}{\int }_{a_i}^{x_i}\frac{f(x_1,...,x_{i-1},u,x_{i+1},...,x_d)}{(x_i-u{)}^{1+\alpha -n}}du$$

$$\frac{{}^C{\partial }^{\alpha }}{\partial x_i^{\alpha }}f(x_1,...,x_d)=\frac{1}{\Gamma (n-\alpha )}{\int }_{a_i}^{x_i}\frac{{\partial }_u^{n}f(x_1,...,x_{i-1},u,x_{i+1},...,x_d)}{(x_i-u{)}^{1+\alpha -n}}du$$

$$u\frac{\partial \bar{c}}{\partial x}=\frac{\partial }{\partial x}(K_x\frac{\partial \bar{c}}{\partial x})+\frac{\partial }{\partial y}(K_y\frac{\partial \bar{c}}{\partial y})+\frac{\partial }{\partial z}(K_z\frac{\partial \bar{c}}{\partial z}),$$

$$u\frac{\partial \overline{c^y}}{\partial x}=\frac{\partial }{\partial z}(K_z\frac{\partial \overline{c^y}}{\partial z}).$$

$$\overline{c^y}(x,z)=\frac{1}{\sqrt{4\pi K_{z}x}}\exp \left(\frac{-z^2}{4K_{z}x}\right).$$

$$u\frac{\partial \overline{c^y}(x,k)}{\partial x}=-K_z{k}^2\overline{c^y}(x,k),$$

$$\overline{c^y}(x,k)=\overline{c^y}(0,k)\exp \left(-K_z{k}^{2}x\right),$$

$$⟨z^2(x)⟩=-\underset{k\to 0}{\lim }\frac{d^2}{d{k}^2}\overline{c^y}(x,k)=2K_{z}x.$$

$$\overline{c^y}(s,k)=\frac{\overline{c^y}(0,k)}{s+K_z{s}^{1-\alpha }{k}^2},$$

$$⟨z^2(x)⟩=-{\mathcal{L}}^{-1}\left\{\underset{k\to 0}{\lim }\frac{d^2}{d{k}^2}\overline{c^y}(s,k)\right\}=\frac{2K_z}{\Gamma (\alpha +1)}{x}^{\alpha },$$

$$\mathcal{L}\left\{{}_0{J}_x^{\alpha }f(x)\right\}=s^{-\alpha }f(s),$$

$$\overline{c^y}(x,z)-\overline{c^y}(0,z)={}_0{J}_x^{\alpha }{K}_z\frac{{\partial }^2\overline{c^y}(x,z)}{\partial z^2},$$

$$\frac{\partial \overline{c^y}(x,z)}{\partial x}=\frac{{\partial }^{1-\alpha }}{\partial x^{1-\alpha }}{K}_z\frac{{\partial }^2\overline{c^y}(x,z)}{\partial z^2}.$$

$$\frac{{\partial }^{\alpha }\overline{c^y}(x,z)}{\partial x^{\alpha }}-\frac{x^{-\alpha }}{\Gamma (1-\alpha )}\overline{c^y}(0,z)=K_z\frac{{\partial }^2\overline{c^y}(x,z)}{\partial z^2}.$$

$${}_0^C{D}_x^{\alpha }f(x)={}_0{D}_x^{\alpha }f(x)-\frac{f(0)}{\Gamma (1-\alpha )}{x}^{-\alpha }$$

$$u\frac{{}^C{\partial }^{\alpha }\overline{c^y}(x,z)}{\partial x^{\alpha }}=K_z\frac{{\partial }^2\overline{c^y}(x,z)}{\partial z^2},$$

$$u(z)=\gamma z^k,\gamma =\frac{u(z_1)}{z_1^k}$$

$$\gamma \frac{{}^C{\partial }^{\alpha }\overline{c^y}}{\partial x^{\alpha }}=K\frac{{\partial }^2\overline{c^y}}{\partial z^2}.$$

$$K_z=K(x)=\rho x^{\alpha },\rho =\nu u,$$

$$\gamma z^k\frac{{}^C{\partial }^{\alpha }\overline{c^y}}{\partial x^{\alpha }}=\rho x^{\alpha }\frac{{\partial }^2\overline{c^y}}{\partial z^2}.$$

$$K_z\frac{\partial \overline{c^y}}{\partial z}=0,z=z_0,z=h,$$

$$u\overline{c^y}(0,z)=Q\delta (z-H_s),x=0,$$

$$\overline{c^y}(x,z)=X(x)Z(z).$$

$${}_0^C{D}_x^{\alpha }X+\kappa {\lambda }^{2}X=0,$$

$$\frac{d^{2}Z}{d{z}^2}+{\lambda }^{2}Z=0.$$

$${}_0^C{D}_x^{\alpha }X+\tau {\lambda }^2{x}^{\alpha }X=0,$$

Full text

I Introduction

The dispersion of pollutants in the atmosphere is a permanent source of challenging problems due to its physical complexity. An example of nontrivial problems is the description of diffusion under atmospheric turbulence. Indeed, the turbulence is the reason behind the dispersion of pollutants in the atmosphere, since without turbulence the pollutants would follow only the streamlines of mean wind velocities, displaying minimal diffusion in other directions. A remarkable consequence of turbulence is the emergence of anomalous diffusion. Unlike common diffusion, where the mean square displacement increases linearly with time, in the anomalous diffusion the mean square displacement is not linear. The anomalous diffusion is closely connected with the failure of the central limit theorem due to sparse distribution or long-range correlations. Actually, the anomalous diffusion is related to the more general Lévy–Gnedenko theorem that generalizes the central limit theorem for situations where not all moments exist Metzler . Historically, anomalous diffusion was first observed in nature in the pollutant dispersion phenomenon. In 1926, Richardson measured the increase in the width of plumes of smoke generated from point sources located in a turbulent velocity field Metzler ; Metzler2 ; Richardson ; West . Based on its observations, Richardson speculated that the speed of turbulent air, which has a non-differentiable structure, can be approximately described by the Weierstrass function. This was motivated partly by the observation that the width of smoke plumes grows with $t^{\alpha}$ ($\alpha\geq 3$), unlike common diffusion where $\alpha=1$. Furthermore, the non-differentiable behavior of the width growth of plumes generated from a point source is directly related to the fractal structure of the turbulent velocity field, where the fluctuation’s size scales are, in many cases, very large compared to the average scale.

In this context, traditional differential equations do not adequately describe the problem of turbulent diffusion because usual derivatives are not well defined in the non-differentiable behavior introduced by turbulence. Consequently, it is expected that classical advection-diffusion equations do not fully explain the anomalous diffusion of pollutants since in this case the parameters of the system usually grow faster than the solutions obtained by classical models West . However, despite this facts, traditional Eulerian and Lagrangian models are the most frequently used in modeling the dispersion of contaminants in the atmosphere GL ; HP ; MRVG ; Wilson . The Eulerian models basically consist of the solution of the advection-diffusion equation and the Lagrangian models are based on the solution of the Langevin equation. In order to deal with the anomalous diffusion, a common procedure used in the literature in order to modify Eulerian models is to assume that the physical structure of turbulent flow and velocity fields are described by a complex diffusion coefficient and mean velocity profile that are both considered as functions of spatial coordinates. These functions are usually chosen in order to fit experimental data or are obtained from Taylor statistical diffusion theory Taylor ; Batchelor ; DACMT ; GMCT .

The dispersion model that we propose in the present work, unlike the models usually found in the literature, does not solve the advection-diffusion equation as expressed traditionally, but modifies the mathematical structure of this equation to represent more realistically the evolution in space of concentration of contaminants dispersed in a turbulent flow. In this sense, we introduce fractional operators in the equation that governs the distribution of contaminants in the atmosphere. The use of fractional calculus in the modeling of turbulent diffusion is justified by the non-differentiable behavior in the problem and by the presence of anomalous diffusion. In the last decades, thousands of works were carried out with the objective of explaining anomalous diffusion. However, few works deal with the analysis of the validity of models based on classical differential equations, and or the use of unusual differential operators, to describe systems displaying non-differential behavior and/or anomalous dynamics. In this subject, we can highlight the use of fractional calculus that emerged as an valuable mathematical tool to model the time evolution of anomalous diffusion since fractional derivatives arouses naturally in anomalous diffusion process (a detailed physical explanation for the emergence of fractional derivatives in anomalous diffusion models can be found in the reviews Metzler ; Metzler2 ; West ). However, the use of fractional derivatives in the study of steady states systems, and in particular the dispersion of pollutants in the atmosphere, is underexplored. In this context, in order to investigate the potential of application of fractional derivatives to modeling dispersion of pollutants, we propose two very simple fractional differential equation models for the spatial distribution of concentration of a non-reactive pollutant in Planetary Boundary Layer (PBL). We solve these models and we compare the solutions with traditional integer order derivative models and with a real experiment.

The paper is organized in the following way. In Section II we review the basic definitions and properties of Riemann-Liouville and Caputo fractional derivatives, that are needed for formulating our fractional differential equation model. The models are presented and analytically solved in Section III. A numerical comparison with our solution against traditional models and experimental data is done in Section IV. Finally, Section V presents the conclusions.

II Mathematical background

There are several definitions of fractional order derivatives. Theses definitions include the Riemann-Liouville, Caputo, Riesz, Weyl, Grunwald-Letnikov, etc. (see OldhamSpanier ; Kilbas ; SKM ; MillerRoss ; Podlubny ; Diethelm for a review of fractional calculus and fractional differential equations). In this section, we review the basic definitions and properties of the Riemann-Liouville and Caputo Fractional Calculus. Despite we have many different approaches to fractional calculus, several known formulations are somehow connected with the analytical continuation of Cauchy formula for $n$-fold integration, that give us a definition for an integration $\!\!\!\!\!\phantom{J}{}_{a}J^{\alpha}_{x}$ of noninteger (or fractional) order defined by:

[TABLE]

with $a<b$ and $a,b\in\mathbb{R}$. This integral is historically called left fractional Riemann-Liouville integral of order $\alpha\in\mathbb{R}$. For integer $\alpha$ the fractional Riemann-Liouville integral (1) coincides with the usual integer order $n$-fold integration. Moreover, from the definition (1) it is easy to see that the Riemann-Liouville fractional integral converges for any integrable function $f$ if $\alpha>1$. Furthermore, it is possible to proof the convergence of (1) for $f\in L_{1}[a,b]$ even when $0<\alpha<1$ Diethelm .

The fractional order integration (1) is the building block of the Riemann-Liouville and Caputo calculus, the two most popular formulations of fractional calculus, as well as several other approaches to fractional calculus. The left Riemann-Liouville fractional derivative of order $\alpha>0$ ($\alpha\in\mathbb{R}$) is defined by ${{}_{a}D^{\alpha}_{x}}f(x)=D^{n}_{x}\!\!\!\!\phantom{J}{}_{a}J^{n-\alpha}_{x}f(x)$ with $n=[\alpha]+1$, namely:

[TABLE]

On the other hand, the left Caputo fractional derivative of order $\alpha>0$ ($\alpha\in\mathbb{R}$) is defined by ${{}_{a}^{C}D^{\alpha}_{x}}f(x)=\!\!\!\!\phantom{J}_{a}J^{n-\alpha}_{x}D^{n}_{x}f(x)$ with $n=[\alpha]+1$, namely:

[TABLE]

where $f^{(n)}(u)=\frac{d^{n}f(u)}{du^{n}}$ are ordinary derivatives of integer order $n$, and $f^{(n)}\in L_{1}[a,b]$.

An important consequence of definitions (2) and (3) is that the Riemann-Liouville and Caputo fractional derivatives are nonlocal operators. The left differ-integration operators (2) and (3) depend on the values of the function at left of $x$, i.e. $a\leq u\leq x$ ($x\leq u\leq b$). Moreover, when $\alpha$ is an integer the Riemann-Liouville and Caputo derivatives reduce to ordinary derivatives of order $\alpha$. It is important to notice the remarkable fact that while the Riemann-Liouville derivative of a constant is not zero, the Caputo derivatives of a constant is zero for all $\alpha>0$. Furthermore, for power functions $(x-a)^{\beta}$ and $(b-x)^{\beta}$ ($\beta>0$) we have:

[TABLE]

Finally, let us consider now a real valued function $f:\Omega\rightarrow\mathbb{R}$ of $d$ real variables $x_{1},x_{1},...x_{d}$ defined over the domain $\Omega=[a_{1},b_{1}]\times\cdots\times[a_{d},b_{d}]\subset\mathbb{R}^{d}$. We can define the left Riemann-Liouville and Caputo partial fractional derivatives of order $\alpha\in\mathbb{R}_{+}^{*}$ with respect to $x_{i}$ ($n=[\alpha]+1$, $a_{i},b_{i}\in\mathbb{R}$) as

[TABLE]

and

[TABLE]

respectively, where $\partial_{u}^{n}$ is the ordinary partial derivative of integer order $n$ with respect to the variable $u$.

III The mathematical models

An equation for the spatial distribution of concentration of a non-reactive pollutant in PBL can be obtained by an application of the principle of continuity or conservation of mass, where the flows are represented by the K-theory Csanady ; Blackadar . In a Cartesian coordinate system in which the longitudinal direction $x$ coincides with the average wind velocity, the spatial distribution of concentration $\bar{c}=\bar{c}(x,y,z)$ of a non-reactive substance can be described by the advection-diffusion equation

[TABLE]

where we consider only the steady state case, $u=u(z)$ is the mean wind speed in the longitudinal direction and $K_{x}$, $K_{y}$, $K_{z}$ are the diffusion coefficient. The equation of the cross-wind integrated concentration ($\overline{c^{y}}=\overline{c^{y}}(x,z)$) is obtained by integrating the Eq. (7) with respect to $y$ from $-\infty$ to $+\infty$ (neglecting the longitudinal diffusion),

[TABLE]

A process governed by the steady state regime of the advection-diffusion equation (8) with $K_{z}$ constant is called a Gaussian process with a sharp boundary condition $\overline{c^{y}}_{0}(z)=\lim_{x\rightarrow 0^{+}}\overline{c^{y}}(x,z)=\delta(z)$ in the domain $0\leq x<\infty$ and $-\infty<z<\infty$ the solution of (8) is given by the Gaussian shape

[TABLE]

On the other hand, by taking a Fourier transform of (8) we obtain a relaxation equation for a fixed wavenumber $k$,

[TABLE]

which implies that individual modes of (8) on rectangular domains decay exponentially in $x$

[TABLE]

where we set for simplicity $u=1$. Consequently, for a Gaussian process, individual modes of the concentration follows an exponential distribution in the longitudinal direction $x$. Actually, the Gaussian distribution for the concentration (9) and its individual modes exponential distribution (11) are characteristic of the normal diffusive process, displaying a linear mean squared displacement

[TABLE]

On the other hand, in most cases, the anomalous process displays a power law mean squared displacement $\langle z^{2}\rangle\propto x^{\alpha}$ with $0<\alpha<1$ Metzler . This power law follows directly from the Fourier-Laplace concentration function Metzler

[TABLE]

since

[TABLE]

where $\Gamma(\cdot)$ is the Gamma function and $\mathscr{L}^{-1}$ stands for the inverse Laplace transform. Moreover, from the Laplace transform for fractional Riemann-Liouville integrals OldhamSpanier ; SKM

[TABLE]

we can infer from (13) the steady-state fractional advection-diffusion integral equation

[TABLE]

for $u=1$. Here and in the rest of the work we set without loss of generality $a=0$. Furthermore, by taking the derivative of (16) in $x$ we obtain a Riemann-Liouville fractional differential equation

[TABLE]

Although it is well established that fractional differential equations describe more realistically the time evolution of anomalous diffusive systems Metzler ; Metzler2 ; West , in the present work we shall show that the fractional differential equation (17) is more adequate than the classical equation (8) to describes the steady state regime of anomalous dispersion of pollutants in a turbulent flow. Before proceeding the work, let us rewrite (17) in a more adequate expression. Given that for continuum functions $f(x)$ we have ${{}_{0}J^{\alpha}_{x}}{{}_{0}D_{x}^{\alpha}}f(x)={{}_{0}D_{x}^{\alpha}}{{}_{0}J^{\alpha}_{x}}f(x)=f(x)$, and ${{}_{0}J^{\alpha}_{x}}{{}_{0}J^{\beta}_{x}}={{}_{0}J^{\alpha+\beta}_{x}}$ OldhamSpanier ; Kilbas ; SKM ; Diethelm , we obtain after some manipulations

[TABLE]

Finally, the Riemann-Liouville derivative in (18) can be replaced by a Caputo derivative since they are directly related

[TABLE]

for $0<\alpha<1$ Kilbas ; Diethelm . Consequently, the general steady-state fractional advection-diffusion equation displaying the anomalous power law mean squared displacement (14) is given by

[TABLE]

for $u$ constant. The practical advantage of (19) instead of the equivalent equation (17) is that while differential equations containing Caputo derivatives requires regular physical boundary conditions, differential equations containing Riemann-Liouville derivatives requires non-regular boundary conditions that are more difficult to implement and to physically interpret HeymansPodlubny .

Actually, in traditional atmospheric dispersion models, the parameters $u$ and $K_{z}$ are usually considered as functions on $x$ and $z$ in order to deal with anomalous diffusion introduced by turbulence. In this context, the turbulent parameterization is related to the estimate of the mean wind speed and the diffusion coefficient Moreira ; Sharan . The wind speed is obtained from the similarity theory Panofsky ; Smith ; Irwin and the diffusion coefficient is obtained from Taylor statistical diffusion theory Taylor ; Batchelor . In this work is used the power law profile for mean wind speed Smith ; Irwin ,

[TABLE]

where $u(z_{1})$ is the wind speed at known height $z_{1}$ and $k$ is the power law exponent dependent on the intensity of turbulence. In order to obtain a simple analytic solution of equation (18), the diffusion coefficient can be specified as a function of downwind distance. We consider here two important cases. The first when $u=\gamma$ and $K_{z}=K$ are constants give us a fractional generalization of the Gaussian model (that we shall call $\alpha$–Gaussian model)

[TABLE]

The second case is a fractional generalization of the model proposed by Sharan and Modani (SM model) Sharan , where now we have

[TABLE]

where $\nu=(\frac{\sigma_{w}}{u})^{2}$ is the turbulence parameter Arya , $\sigma_{w}$ is the standard deviation of the vertical component of the wind speed, and $0<\alpha<1$ is the same as the order of the fractional derivative in (8). Substituting (20) and (22) into (18) we define the $\alpha$–SM model

[TABLE]

It is important to remark that the diffusion coefficient $K_{z}=K(x)=\rho x$ considered in the SM model Sharan (and (22) in $\alpha$–SM model) is not physically meaningful in the context of diffusion of pollutants in real atmospheric condition. This diffusion coefficient is introduced in Sharan only in order to fit experimental data with a simple classical advection-diffusion equation. In the present work, we analyzed the $\alpha$–SM model (23) with (22) only in order to investigate if there is any advantage, against the $\alpha$–Gaussian model (21), of considering the diffusion coefficient as a function of spatial coordinates.

Finally, in order to (18), (21) and (23) describe a possible real problem of dispersion in PBL, it should be imposed boundary conditions of zero flux on the ground ($z=z_{0}$) and top ($z=h)$, and the pollutant is released from an elevated point source of emission rate $Q$ at height $H_{s}$, i.e.,

[TABLE]

[TABLE]

where $z_{0}$ is the surface roughness length, $h$ is the PBL height and $\delta(\cdot)$ is the Dirac delta function.

Equations (21) and (23) can be analytically solved by using the separation of variables technique. We begin by assuming that both (21) and(23) have a solution of the form

[TABLE]

Substituting (26) in (21), we get two ordinary differential equations in the variables $X$ and $Z$ as follows,

[TABLE]

and

[TABLE]

where $\kappa=\frac{K}{\gamma}$. On the other hand, by substituting (26) in (23) we get

[TABLE]

and

[TABLE]

where $\tau=\frac{\rho}{\gamma}$.

III.1 The solution of the $\alpha$–Gaussian model (21)

The solution of (21) subject to the boundary conditions (24) and (25) can be obtained from (27) by the Frobenius method. The solution for the $\alpha$–Gaussian model (21) with $z_{0}\approx 0$ is (see Appendix A)

[TABLE]

where $\lambda_{n}=\frac{n\pi}{h}$, and

[TABLE]

is the Mittag-Leffler function. The solution for the classical Gaussian model is obtained from equation (21) when $\alpha=1$ Csanady

[TABLE]

The big advantage of fractional differential equations, against integer order ones, to describe anomalous diffusion can easily be seen from the solution of the $\alpha$-Gaussian model (31). Even for a constant diffusion coefficient, the concentration of pollutants in (31) falls more slowly in $x$ compared with the traditional Gaussian model (33), since the Mittag-Leffler function in (31) for $0<\alpha<1$ falls more slowly than the exponential function in (33). Actually, this was already expected since the $\alpha$–Gaussian model displays a power law mean squared displacement (14). Furthermore, the function $E_{\alpha}(-\kappa\lambda_{n}^{2}x^{\alpha})$ appearing in the solution of the $\alpha$-Gaussian model (31) is just related to the Mittag-Leffler distribution $1-E_{\alpha}(-\kappa\lambda_{n}^{2}x^{\alpha})$ that is characteristic of anomalous diffusion process like the fractional Brownian motion Metzler ; Metzler2 (it is also well known that the Lévy distribution, found in Lévy flight anomalous diffusion, and Mittag-Leffler distribution are related Feller ; HMS ).

III.2 The solution of the $\alpha$–SM model (23)

The solution of the fractional $\alpha$–SM model (23) subject to the boundary conditions (24) and (25) can be obtained from (29) and (30) by the Frobenius method. The solution of (29) is (see Appendix B)

[TABLE]

On the other hand, the solution of (30) is (Appendix B),

[TABLE]

where $J_{\mu}$ and $J_{-\mu}$ are the Bessel functions of the first kind and order $\mu$ and $-\mu$, respectively, and $\mu=\frac{1}{2+k}$. Applying the boundary condition given by (24) at $z=z_{0}$ (with $z_{0}\approx 0$), we get (Appendix B)

[TABLE]

From the boundary condition given by (24) at $z=h$, we obtain the following relationship (Appendix B),

[TABLE]

Therefore, the $\lambda=\lambda_{n}$ are the zeros of the Bessel function of order $-\mu+1$. Substituting (34), (35), (36) and (37) in (26) we obtain the equation for the cross-wind integrated concentration of a non-reactive pollutant in PBL released from an elevated point source

[TABLE]

The constants $B_{0}$ and $B_{n}$ are obtained from the boundary condition given by (25), and $\lambda_{n}$ are given by the zeros of the Bessel function of order $-\mu+1$ calculated from (37). The function (38) is the solution of our generalized model (23) obtained when $u$ and $K_{z}$ are functions of the $z$ ( (20) and (22)).

IV Performance of the models

The performance of the fractional derivative models expressed in (31) and (38) was evaluated against experimental ground-level concentration using the Copenhagen dispersion experiment Briggs1973 . The experiment took place in Copenhagen in $1978$. The tracer was released without buoyancy at a height of $115m$ and collected at ground level positions at a maximum of three crosswind arcs of tracer sampling units. The site was mainly residential with a roughness length of the $0.6m$ Gryning ; Gry-Lyck . The meteorological conditions of the Copenhagen dispersion experiment are shown in Table $1$, where $\bar{u}_{10}$ is the mean wind speed at $10m$, $u_{*}$ is the friction velocity and $L$ is the Monin–-Obukhov length.

Usually, the performance of dispersion models is evaluated from a well know set of statistical indices described by Hanna Hanna defined in the following way,

[TABLE]

where $c_{p}$ is the computed concentration, $c_{o}$ is the observed concentration, $\sigma_{p}$ is the computed standard deviation, $\sigma_{o}$ is the observed standard deviation, and the overbar indicates an averaged value. The statistical index FA2 represents the fraction of data for $0.5\leq\frac{c_{p}}{c_{o}}\leq 2$. The best results are indicated by values nearest to [math] in NMSE, FS and FB, and nearest to $1$ in Cor and FA2. The statistical indices used to evaluate our model’s performance is shown in Tables $3$.

In Table $2$, the results of the estimated crosswind-integrated concentrations for the Copenhagen experiment Briggs1973 obtained for the fractional models are compared with experimental data, with the Gaussian model (33), and with results obtained from the SM model Sharan (where it was proposed a model given by (23) with $\alpha=1$). In Gaussian model (33) and $\alpha$–Gaussian model (31) the mean wind speed in longitudinal direction ($u$) is obtained from Table 1 and the diffusion coefficient constant $K$ is chosen, in order of comparison, as being the mean of the diffusion coefficient given by equation (22) (with $\alpha=1$) over the longest distance $x$ of each experiment.

The good performance of the proposed fractional models is shown in Tables $2$, $3$. In special, Table $3$ shows the statistical indices obtained from the mathematical models and the results obtained in the Copenhagen dispersion experiment Briggs1973 . The $\alpha$–Gaussian model (31) performs far better than the Gaussian model (33) to describe the data of Copenhagen experiment, as can be seen from the statistical indices. The advantage of the $\alpha$–Gaussian model is a consequence of the anomalous diffusion present on turbulent diffusion that results in a distribution with a power law mean squared displacement. It is important to remark, that the $\alpha$–Gaussian model also has better performance than the SM model Sharan , where it was considered that the diffusion coefficient is a function $K_{z}=K(x)=\rho x$ in order to deal with the anomalous diffusion of the problem. Furthermore, from Table 3 we see that the $\alpha$–SM model (38) has a very small extra performance over the $\alpha$–Gaussian model. This result indicates that there may not be much advantage to consider the diffusion coefficient as a not physically meaningful function. Finally, in order to estimate the better $\alpha$ value for each model, we analyzed the solutions from $\alpha=0.7$ to $\alpha=0.99$ by steps of $0.01$. The results obtained from the $\alpha$–SM (38) and shown in Table 3 are for $\alpha=0.92$. For this value of $\alpha$, the model (38) has its best statistical indices. For the $\alpha$–Gaussian model (31), the best performance occurs when $\alpha=0.80$. This difference in the best $\alpha$–value is related to different diffusion coefficients employed in both cases. In the Gaussian model (33) the diffusion coefficient ($K$) is constant, but in $\alpha$–SM model the diffusion coefficient is a function of the distance $x$ (38). This suggests a correlation between the order of the equation and the quality of diffusion coefficient to the model best describe the observed concentration data. The traditional models ($\alpha=1$) originate from the equation of molecular diffusion (Fick’s law), which assumes a Gaussian distribution that displays a linear mean square displacement. Asymmetries related to turbulent flow in the atmosphere are described by the diffusion coefficient. Already the fractional model assumes an anomalous probability distribution with a power law mean squared displacement (14). This anomalous distribution has been proved to be more efficient to describe the motion of particles in a turbulent flow Metzler .

V Conclusion

The main objective of the present work is to investigate the potential application of fractional derivatives to model diffusion of pollutants in the atmosphere. The use of fractional calculus in modeling diffusion of pollutants is justified by the presence of anomalous diffusion due to turbulence. In the last decades, thousands of works were carried out with the objective of explaining anomalous diffusion. However, in this subject only a few works deal with the analysis of the validity of models based on classical differential equations, and or the use of unusual differential operators, to describe systems displaying non-differential behavior and/or anomalous dynamics. On the other hand, in the last decades, the use of fractional differential operators emerged as a valuable mathematical tool to model the time evolution of anomalous diffusion Metzler ; Metzler2 ; West . However, the use of fractional derivatives to study steady states regimes of dispersion of pollutants in the atmosphere is actually underexplored. In this context, in order to investigate the potential of application of fractional derivatives to model dispersion of pollutants, we propose two simple fractional differential equation models for the spatial distribution of concentration of a non-reactive pollutant in PBL. These equations are obtained by considering a system displaying a distribution with a power law mean squared displacement, instead of a linear one as in normal diffusion. We solve these models and we compare the solutions with traditional integer order derivative models and with a real experiment. We found that the fractional derivative $\alpha$–Gaussian model (31) performs far better than the Gaussian model and even better than the SM model Sharan , where it was considered that the diffusion coefficient is a function of the position in order to deal with the anomalous diffusion. Most important, the comparison between the $\alpha$–Gaussian model (31) and the fractional derivative $\alpha$–SM model (38) indicates that, by using fractional derivatives, there can be little advantage to consider the diffusion coefficient as a function of position, since the $\alpha$–Gaussian model with constant diffusion coefficient gives a very good performance.

Finally, the results obtained from the fractional derivative models motivates further investigations of applications of fractional differential equations to model diffusion of pollutants.

Acknowledgments

This work was supported in part by CNPq and CAPES, Brazilian funding agencies.

Appendix A Solution of the $\alpha$–Gaussian model

In order to solve the equation (27) by Frobenius method, we consider the following general power series

[TABLE]

where $0<\alpha<1$, and $p,\beta$ are constants. The Caputo fractional derivative of (39) is then given by,

[TABLE]

Without loss of generality, by considering $\beta=p+1$ we have,

[TABLE]

where,

[TABLE]

If $p<0$, the function $\frac{dx^{p}}{dx}$, and consequently $\frac{dX(x^{\alpha})}{dx}$, is not integrable and the Caputo fractional derivative is not defined. Substituting equation (41) into equation (27) we obtain

[TABLE]

If $p\neq 0$ we have the trivial solution $X(x^{\alpha})=0$ since for $b_{0}\neq 0$ the function is not integrable (because in this case $\beta=0,-1,-2,...$ and $p=\beta-1<0$ in order to satisfy (43)). For $p=0$, we obtain a nontrivial solution

[TABLE]

where $E_{\alpha}$ is the Mittag-Leffler function (32).

On the other hand, the solutions $Z_{n}(z)$ ($n=0,1,2,3,...$) of (28) that satisfy the boundary conditions (24) with $z_{0}=0$ are $Z_{n}(z)=B_{n}\cos(\lambda_{n}z)$ where $\lambda_{n}=\frac{n\pi}{h}$, and $B_{n}$ is a constant. Consequently, starting from the superposition principle, and from (26), we obtain

[TABLE]

Finally, by imposing the boundary condition (25), and using the identity

[TABLE]

we obtain (31).

Appendix B Solution of the $\alpha$–SM model

In order to solve the equations (29) and (30), we use Frobenius method. Substituting equation (41) in equation (29) we obtain,

[TABLE]

As in Appendix A, if $p\neq 0$ the Caputo fractional derivative for the function $X(x^{\alpha})$ is not defined. Therefore, $p=0$ and $\beta=1$. In this case, the equation (47) can be written as,

[TABLE]

From the equation (48) we have $b_{1}=0$ and $b_{2j+1}=0$. Furthermore, for $b_{2j}$ we get

[TABLE]

By substituting (49) in (39) (for $b_{0}=1)$ we get the equation (34),

[TABLE]

On the other hand, (30) can also be solved by Frobenius method. We consider the following power series as analytic solution of (30)

[TABLE]

Substituting $Z(z)$ into (30) we have (shifting the indices of first two series so all terms are of the form $z^{n+r}$),

[TABLE]

All coefficients of powers \Big{(}z^{1+\frac{k}{2}}\Big{)}^{n+r} must equate to zero to obtain a solution. This yields two cases. In the first case ($c_{1}=0$ and $r=\displaystyle\frac{2}{2+k}$) the solution is given by,

[TABLE]

In the second case ($c_{0}=0$ and $r=-1$), the solution is

[TABLE]

Using the superposition principle we have that

[TABLE]

where $A$ and $B$ are constants. Applying the boundary condition $Z^{\prime}(0)=0$ (24) in (55) we obtain $A=0$. On the other hand, by applying the boundary condition $Z^{\prime}(h)=0$ (24) in (55) with $A=0$ we obtain

[TABLE]

Making the change of variable $m=n-1$ we get

[TABLE]

Therefore we should have

[TABLE]

Consequently, $\lambda=\lambda_{n}$ is given by the zeros of the Bessel function of order $-\mu+1$.

Finally, from the superposition principle, and from (26), (50) and (55), the solution of (23) can be written as

[TABLE]

In order to fix the constants $B_{n}$, we use the boundary condition $C(0,z)=\frac{Q}{\gamma z^{k}}\delta(z-H_{s})$ (25) in (59). We consider

[TABLE]

as a good approximation because $J_{-\mu}(2\lambda_{n}\mu z^{\frac{1}{2\mu}})$ goes to zero when $n\rightarrow\infty$. Now, making

[TABLE]

where $l\in\{0,1,...,m\}$, we fix $B_{n}$ from the numerical solution of the following system

[TABLE]