A new parameterization for the concentration flux using the fractional calculus to model the dispersion of contaminants in the Planetary Boundary Layer

TL;DR

This paper introduces a fractional calculus-based parameterization for modeling contaminant dispersion in the Planetary Boundary Layer, demonstrating superior fit to experimental data over traditional models.

Contribution

It develops a fractional differential equation model for concentration flux that outperforms classical models and suggests a link between fractional order and turbulence structure.

Findings

Fractional model fits experimental data better than Gaussian models.

Constant parameter fractional model outperforms some variable-parameter models.

The fractional order relates to the physical turbulence structure.

Abstract

In the present work, we propose a new parameterization for the concentration flux using fractional derivatives. The fractional order differential equation in the longitudinal and vertical directions is used to obtain the concentration distribution of contaminants in the Planetary Boundary Layer. We solve this model and we compare the solution against both real experiments and traditional integer order derivative models. We show that our fractional model gives very good results in fitting the experimental data, and perform far better than the traditional Gaussian model. In fact, the fractional model, with constant wind speed and a constant eddy diffusivity, performs even better than some models found in the literature where it is considered that the wind speed and eddy diffusivity are functions of the position. The results obtained show that the structure of the fractional order differential equation is more appropriate to calculate the distribution of dispersed contaminants in a turbulent flow than an integer-order differential equation. Furthermore, a very important result we found it is that there should be a relation between the order $\alpha$ of the fractional derivative with the physical structure of the turbulent flow.

Tables

Table 1. Table 1: Copenhagen Experiment for h | L | < 10 ℎ 𝐿 10 \frac{h}{|L|}<10 (instable)

Model	Cor	NMSE	FS	FB	FA2
Eq.(11)	0.97	0.05	0.08	-0.24	1.00
Gauss	0.96	0.17	0.06	-0.44	0.75
O-G	0.97	0.83	1.00	-0.77	0.41
Moreira (2005)	0.97	0.02	0.05	0.01	1.00
Kumar (2012)	0.90	0.05	0.34	-0.04	0.96

Table 2. Table 2: Prairie Grass Experiment for h | L | < 4 ℎ 𝐿 4 \frac{h}{|L|}<4 (stable)

Distance (m)	Model	Cor	NMSE	FS	FB	Fa2
50	Eq.(11)	0.96	4.39	1.22	-1.47	0.00
	Gauss	0.96	4.91	1.53	-1.47	0.00
	O-G	0.96	4.91	1.53	-1.47	0.00
200	Eq.(11)	0.96	1.68	1.06	-1.08	0.00
	Gauss	0.96	0.95	1.08	-0.83	0.06
	O-G	0.96	0.95	1.08	-0.83	0.06
800	Eq.(11)	0.90	0.47	0.89	-0.56	0.88
	Gauss	0.84	1.72	1.35	-1.02	0.06
	O-G	0.84	1.72	1.35	-1.02	0.06
$\ge 200$	Eq.(11)	0.91	1.35	1.20	-0.89	0.41
	Gauss	0.97	2.82	1.46	-1.20	0.06
	O-G	0.97	2.82	1.46	-1.20	0.06
$\ge 400$	Eq.(11)	0.93	0.73	1.03	-0.68	0.41
	Gauss	0.95	2.06	1.40	-1.08	0.08
	O-G	0.95	2.06	1.40	-1.08	0.08
All distances	Eq.(11)	0.85	3.72	1.47	-1.27	0.31
	Gauss	0.98	4.56	1.54	-1.39	0.02
	O-G	0.98	4.56	1.54	-1.39	0.02

Table 3. Table 3: Hanford Experiment for h | L | ≤ 3 ℎ 𝐿 3 \frac{h}{|L|}\leq 3 (stable)

Distance (m)	Model	Cor	NMSE	FS	FB	Fa2
100	Eq.(11)	0.65	10.32	1.80	-1.66	0.00
	Gauss	0.63	13.21	1.84	-1.71	0.00
	O-G	0.62	15.33	1.88	-1.74	0.00
1600	Eq.(11)	0.97	0.36	0.72	-0.53	1.00
	Gauss	0.97	0.42	0.69	-0.58	0.66
	O-G	0.88	6.21	1.78	-1.49	0.00
3200	Eq.(11)	0.95	0.10	0.50	-0.19	1.00
	Gauss	0.95	0.12	0.49	-0.22	1.00
	O-G	0.84	6.54	1.81	-1.48	0.00
$\le 800$	Eq.(11)	0.60	8.27	1.77	-1.51	0.05
	Gauss	0.48	10.32	1.79	-1.58	0.00
	O-G	0.82	14.24	1.83	-1.71	0.00
$\ge 1600$	Eq.(11)	0.92	0.26	0.67	-0.38	1.00
	Gauss	0.90	0.29	0.65	-0.42	0.83
	O-G	0.83	6.31	1.75	-1.49	0.00
All distances	Eq.(11)	0.55	7.57	1.74	-1.35	0.43
	Gauss	0.39	9.31	1.76	-1.42	0.33
	O-G	0.87	14.01	1.80	-1.68	0.00

Table 4. Table 4: Copenhagen Experiment for h | L | > 10 ℎ 𝐿 10 \frac{h}{|L|}>10 (instable)

Model	Cor	NMSE	FS	FB	FA2
Eq.(11)	0.65	0.20	0.97	-0.14	0.90
Gauss	0.62	0.34	1.02	-0.33	0.81
O-G	0.71	0.47	1.09	-0.47	0.72

Table 5. Table 5: Prairie Grass Experiment for h | L | > 10 ℎ 𝐿 10 \frac{h}{|L|}>10 (instable)

Distance (m)	Model	Cor	NMSE	FS	FB	Fa2
50	Eq.(11)	0.80	0.37	-1.04	0.07	0.75
	Gauss	0.89	0.40	-0.76	0.58	0.75
	O-G	0.89	0.40	-0.76	0.58	0.75
100	Eq.(11)	0.48	1.26	-1.59	0.66	0.65
	Gauss	0.41	1.24	-1.38	0.91	0.20
	O-G	0.41	1.24	-1.38	0.91	0.20

Equations

$$\frac{D\bar{c}}{Dt}+\overrightarrow{\nabla }\cdot {\overrightarrow{\Pi }}_c=0,$$

$${\overrightarrow{\Pi }}_c=-K\overrightarrow{\nabla }\bar{c},$$

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

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

$${\Pi }_{c,z}=-K_z\frac{{\partial }^{\beta }\overline{c}}{\partial z^{\beta }},$$

$$\frac{{\partial }^{\beta }f(x,y,z)}{\partial z^{\beta }}=\frac{1}{\Gamma (1-\beta )}{\int }_a^z\frac{{\partial }_{u}f(x,y,u)}{(z-u{)}^{\beta }}du,$$

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

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

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

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

$$\overline{c^y}(x,z)=\sum _{n=0}^{\infty }{a}_n{E}_{\alpha }(-\kappa {\lambda }_n^2{x}^{\alpha })E_{\alpha +1}(-{\lambda }_n^2{z}^{\alpha +1}),$$

$$E_{\alpha }(x)=\sum _{n=0}^{\infty }\frac{x^n}{\Gamma (n\alpha +1)},$$

$$\overline{c^{y}}(x,z)=\frac{Q}{uh}\Big{[}1+2\sum_{n=1}^{\infty}\cos(\lambda_{n}H_{s})\cos(\lambda_{n}z)\exp(-\kappa\lambda_{n}^{2}x)\Big{]}.$$

$$\overline{c^{y}}(x,z)=\frac{Q}{2\sqrt{\pi\kappa x}u}\Big{[}\exp\Big{(}-\frac{(z-H_{s})^{2}}{4\kappa x}\Big{)}+\exp\Big{(}-\frac{(z+H_{s})^{2}}{4\kappa x}\Big{)}\Big{]}.$$

$$\begin{array}{rl}\text{NMSE}(\text{normalized mean square error})& =\frac{\overline{(c_o-c_p{)}^2}}{\overline{c_o}\overline{c_p}},\\ \text{Cor}(\text{correlation coefficient})& =\frac{\overline{(c_o-\overline{c_p})(c_p-\overline{c_p})}}{{\sigma }_o{\sigma }_p},\\ \text{FB}(\text{fractional bias})& =\frac{\overline{c_o-\overline{c_p}}}{0.5(\overline{c_o}+\overline{c_p})},\\ \text{FS}(\text{fractional standard deviations})& =\frac{{\sigma }_o-{\sigma }_p}{0.5({\sigma }_o+{\sigma }_p)},\end{array}$$

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

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

$$\frac{d}{dz}[{}_0^C{D}_z^{\alpha }Z]+{\lambda }^{2}Z=0.$$

$$\begin{array}[]{lll}Z(z^{\alpha})=z^{p}\displaystyle\sum_{j=0}^{\infty}b_{j}\frac{z^{j(\alpha+1)}}{\Gamma((\alpha+1)j+\beta)},\end{array}$$

$$\begin{array}[]{lll}{{}_{0}^{C}D^{\alpha}_{z}}Z(z^{\alpha})&=&\frac{b_{0}}{\Gamma(\beta)}{{}_{0}^{C}D^{\alpha}_{z}}z^{p}+\displaystyle\sum_{j=1}^{\infty}\frac{b_{j}}{\Gamma((\alpha+1)j+\beta)}\frac{\Gamma((\alpha+1)j+p+1)}{\Gamma((\alpha+1)j+p-\alpha+1)}z^{(j-1)\alpha+j+p}.\end{array}$$

$$\begin{array}[]{lll}{{}_{0}^{C}D^{\alpha}_{z}}Z(z^{\alpha})=\frac{b_{0}}{\Gamma(p+1)}D_{z}^{\alpha}z^{p}+\displaystyle\sum_{j=1}^{\infty}\frac{b_{j}}{\Gamma((\alpha+1)j+p-\alpha+1)}z^{(j-1)\alpha+j+p},\end{array}$$

$${{}_{0}^{C}D^{\alpha}_{z}}z^{p}=\left\{\begin{array}[]{cl}&0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;p=0\\ &\frac{\Gamma(p+1)}{\Gamma(p+1-\alpha)}z^{p-\alpha},\;\;\;\;\;\;\;\;0<p\leq 1\end{array}\right..$$

$$\begin{array}[]{lll}\frac{d}{dz}{{}_{0}^{C}D^{\alpha}_{z}}Z(z^{\alpha})=\frac{b_{0}}{\Gamma(p+1)}\frac{d}{dz}{{}_{0}^{C}D^{\alpha}_{z}}z^{p}+\displaystyle\sum_{j=0}^{\infty}\frac{b_{j+1}}{\Gamma((\alpha+1)j+1+p)}z^{j\alpha+j+p},\end{array}$$

$$\begin{array}[]{lll}\frac{d}{dz}{{}_{0}^{C}D^{\alpha}_{z}}Z(z^{\alpha})=\frac{b_{0}(p-\alpha)}{\Gamma(p+1-\alpha)}z^{p-1-\alpha}+\displaystyle\sum_{j=0}^{\infty}\frac{b_{j+1}}{\Gamma((\alpha+1)j+1+p)}z^{(\alpha+1)j+p}.\end{array}$$

$$\begin{array}[]{lll}\frac{b_{0}(p-\alpha)}{\Gamma(p+1-\alpha)}z^{p-1-\alpha}+\displaystyle\sum_{j=0}^{\infty}\frac{1}{\Gamma((\alpha+1)j+1+p)}\Big{(}b_{j+1}+\lambda^{2}b_{j}\Big{)}z^{(\alpha+1)j+p}=0,\end{array}$$

$$Z_1(z)=b_0{z}^{\alpha }\sum _{j=0}^{\infty }\frac{(-{\lambda }^2{)}^j{z}^{(\alpha +1)j}}{\Gamma ((\alpha +1)j+\alpha +1)}=b_0{z}^{\alpha }{E}_{\alpha +1,\alpha +1}(-{\lambda }^2{z}^{\alpha +1}).$$

$$\begin{array}[]{lll}{{}_{0}^{C}D^{\alpha}_{z}}Z_{1}(z)=\frac{b_{0}}{\Gamma(1)}+b_{0}z^{\alpha}\displaystyle\sum_{j=1}^{\infty}\frac{(-\lambda^{2})^{j}}{\Gamma(j(\alpha+1)+1)}z^{(\alpha+1)j-\alpha},\end{array}$$

$$\begin{array}[]{lll}{{}_{0}^{C}D^{\alpha}_{z}}Z_{1}(0)=\frac{b_{0}}{\Gamma(1)}\neq 0.\end{array}$$

$$\begin{array}[]{lll}\frac{d}{dz}{{}_{0}^{C}D^{\alpha}_{z}}Z(z^{\alpha})=\displaystyle\sum_{j=0}^{\infty}\frac{b_{j+1}}{\Gamma((\alpha+1)j+1)}z^{j\alpha+j}.\end{array}$$

$$\displaystyle\sum_{j=0}^{\infty}\frac{1}{\Gamma((\alpha+1)j+1)}\Big{(}b_{j+1}+\lambda^{2}b_{j}\Big{)}z^{j\alpha+j}=0,$$

Full text

††journal: Physica A: Statistical Mechanichs and its Applications

1 Introduction

Atmospheric turbulence plays a fundamental role in the dispersion of contaminants in the Planetary Boundary Layer (PBL). In fact, without turbulence, the contaminants would follow only the streamlines of mean wind velocities displaying minimal diffusion in other directions. From the theoretical point of view, turbulence is a permanent source of challenging problems due to its physical and mathematical complexity. As an example of this complexity, a remarkable consequence of turbulence is the emergence of the anomalous diffusion phenomenon. 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, it is related to the more general Lévy-Gnedenko theorem that generalizes the central limit theorem for situations where not all moments exist [1]. Historically, anomalous diffusion was first observed in nature in the dispersion of contaminants phenomenon in 1926, when Richardson measured the increase in the width of plumes of smoke generated from point sources located in a turbulent velocity field [1, 2, 3, 4, 5]. 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$.

The non-differentiable behavior for 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 [8]. In this context, the classical advection-diffusion equation does not fully explain the diffusion of contaminants since the parameters of the system usually grow faster than the solutions obtained by classical models [4, 5]. Moreover, it is expected that 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. Notwithstanding this fact, traditional Eulerian and Lagrangian models are the most frequently used in modeling the dispersion of contaminants in the atmosphere [9, 10, 11, 12]. The Eulerian models consist basically of the advection-diffusion equation and the Lagrangian models are based on 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 eddy diffusivity and mean velocity profile that are both considered as functions of spatial coordinates. These functions or are usually chosen in order to fit experimental data or they are obtained from the Taylor statistical diffusion theory [13, 14, 15, 16].

Despite the difficulties aroused by the non-differentiable behavior introduced by turbulence, in the last decades, few works dealt with the validity of models based on classical differential equations 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 a valuable mathematical tool to model the time evolution of anomalous diffusion (a detailed physical explanation for the emergence of fractional derivatives in anomalous diffusion models can be found in the reviews [1, 2, 4, 5]. However, the use of fractional derivatives to study steady-state systems, and in particular the dispersion of contaminants in the atmosphere, is largely underexplored. In this context, we demonstrate in [17] that the steady-state regime for a spatial concentration distribution of a non-reactive contaminant, in an anomalous diffusion process displaying a power law mean square displacement, is given naturally by a differential equation containing fractional derivatives. Moreover, we show that the fractional model, even in the simplest case of a constant eddy diffusivity, gives very good results when compared with real data and models usually found in the literature with complex eddy diffusivity and mean velocity profile.

The promising result found in our previous work [17] motivates a deeper analysis of the relationship between the turbulence structure of the turbulent flow in the atmosphere and the anomalous diffusion expressed mathematically by a fractional diffusion equation. In the present work, our main objective is to make this deeper analysis by comparing the solution obtained from our fractional diffusion model with the Copenhagen experiment [18], the Prairie Grass experiment [19] and the Hanford experiment [20]. Furthermore, we also generalize the fractional differential model we introduced in [17] by considering a new parameterization for the concentration flux using the fractional calculus. In the work [17], the parameterization for the concentration flux is the usual one, where it is assumed that the turbulence causes a net movement of material in the direction of the concentration gradient at a rate which is proportional to the magnitude of the gradient [21]. However, the emergence of a fractional operator in the advective term due to the anomalous diffusion process [17] leads to a model with inconsistent space dimensions in its different terms. There are two simple ways to correct this inconsistency. The first is by the introduction of an additional scale parameter [1, 2, 4, 5, 6, 7], and the second one is by the introduction of a fractional derivative in the diffusive term by an appropriate parameterization for the concentration flux. In the present work, we adopted the second approach by considering a parameterization for the concentration flux given by a fractional derivative (gradient) of the concentration. The physical justification for this approach is that the turbulent flow causes a net movement of material that reflects the fractal structure of the turbulent vortices, which can be related to the fractional order of the derivative present in the gradient of the flux [5, 22], motivating the use of this fractional parametrization to investigate the relationship between the turbulence structure of the turbulent flow in the atmosphere and the anomalous diffusion. We show that our fractional model with constant velocity field and a constant eddy diffusivity gives very good results when fitting the experimental data, and also when compared to more complex models with integer order derivatives. Most important, a result we found is that there should be a relation between the order $\alpha$ of the fractional derivative and the physical structure of the turbulent flow. Since, regardless of the experiment, when we have the predominance of mechanical energy input in the turbulent flow all experiments were best described with $\alpha=0.72$, and when we have the predominance of energy input by thermal convection the experimental data is better described by $\alpha=0.80$.

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 fractional model is presented and analytically solved in Section III. A numerical comparison of our solution against traditional models and experimental data is done in Section IV. Finally, Section V presents the conclusions.

2 Model description

An equation for the spatial distribution of concentration of a non-reactive contaminant in the PBL can be obtained by an application of the principle of continuity or conservation of mass [21, 23, 24]:

[TABLE]

where $\frac{D}{Dt}=\frac{\partial}{\partial t}+\vec{u}\cdot\overrightarrow{\nabla}$ is the Lagrangian derivative, $\bar{c}=\bar{c}(x,y,z,t)$ is the average concentration, $\vec{u}$ is the wind velocity, and $\overrightarrow{\Pi}_{c}$ is the concentration flux. A traditional closure for concentration flux problem (1) is based on the gradient-transfer approach, which assumes that turbulence causes a net movement of material in the direction of the concentration gradient, at a rate which is proportional to the magnitude of the gradient [21]:

[TABLE]

where $K$ is the eddy diffusivity. By replacing (2) in (1), considering the steady-state case and choosing a Cartesian coordinate system in which the longitudinal direction $x$ coincides with the average wind velocity, and by neglecting the longitudinal diffusion, we get

[TABLE]

Finally, in order to compare the model with experimental data found in the literature, the equation for the cross-wind integrated concentration ($\overline{c^{y}}=\overline{c^{y}}(x,z)$) is obtained by integrating the Eq. (3) with respect to $y$ from $-\infty$ to $+\infty$:

[TABLE]

The Eq. (4) represents a classical model used to estimate the concentration distribution of a contaminant in the Earth’s atmosphere, obtained by considering closures for concentration flux based on the gradient-transfer approach (2).

In the work [17], the parameterization for the concentration flux is the usual one, given by (2), where it is assumed that the turbulence causes a net movement of material in the direction of the concentration gradient at a rate which is proportional to the magnitude of the gradient [21]. However, the emergence of a fractional operator in the advective term due to the anomalous diffusion process [17] leads to a model with inconsistent space dimensions in its different terms. There are two simple ways to correct this inconsistency. The first is by the introduction of an additional scale parameter [1, 2, 4, 5, 6, 7], and the second one is by the introduction of a fractional derivative in the diffusive term by an appropriate parameterization for the concentration flux. In this work, we adopted the second approach by considering a parameterization for the concentration flux given by a fractional derivative of the concentration,

[TABLE]

where [25, 17],

[TABLE]

is the left Caputo partial fractional derivatives of order $\beta$ ($0<\beta<1$) with respect to $z$, $a\in\mathbb{R}$ and $\partial_{u}$ is the ordinary partial derivative of integer order with respect to the variable $u$ . The physical justification for this approach is that the turbulent flow causes a net movement of material that reflects the fractal structure of the turbulent vortices, which can be related to the fractional order of the derivative present in the gradient of the flux [5, 22]. There are many works in the literature that suggest fractal models in order to describe turbulence (for examples, see [34, 35, 36]). In this context, fractional derivatives are in general more adequated to describe the problem than classical derivatives.

On the other hand, we show in [17] that if the process displays anomalous diffusion with a power law mean squared displacement given by $\langle z^{2}\rangle\propto x^{\alpha}$, with $0<\alpha<1$, then the order of the derivative in the advective term ($x$ direction) would be $\alpha$. In this case, the Eq. (4) becomes,

[TABLE]

where the left Caputo partial fractional derivatives of order $\alpha$ ($0<\alpha<1$) with respect to $x$ is defined as in equation (6). Note that for $\alpha=\beta=1$ (7) reduces to the classical integer order diffusion equation (4), since (7) follows from a direct generalization of (4) for anomalous process [17].

The fractional diffusion equation (7) generalizes the model we introduced in [17] by the inclusion of a Caputo fractional derivative in $z$, as a consequence of the new closure for concentration flux introduced in Eq. (5). A further advantage of our present Eq. (7) against the model in [17] is that for $\alpha=\beta$ it has the same dimension on both sides of (7). Considering that (7) is obtained from the law of conservation of mass, it is physically more correct to have the same dimensions on both sides of equality. In this work, we will consider the simplest case for (7), when $u$ and $K_{z}$ are constants. In this case, the Eq. (7) becomes,

[TABLE]

where $\kappa=\frac{K_{z}}{u}$, and we consider also only the particular case where $\alpha=\beta$. In the particular case where $\alpha=1$, Eq. (8) reduces to the well know advection-diffusion equation (Gaussian model) that describes a Gaussian process.

In order to (8) describes a possible real dispersion process in PBL, it should be imposed boundary conditions of zero flux on the ground ($z=z_{0}$) and top ($z=h$), and consider that the contaminant is released from an elevated point source with emission rate $Q$ at height $H_{s}$, i.e.,

[TABLE]

[TABLE]

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

The solution of the fractional differential equation (8), subjected to the boundary conditions (9) and (10), can be analytically obtained. We have (see A)

[TABLE]

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

[TABLE]

and the constants $a_{n}$ and $\lambda_{n}$ are obtained from the boundary conditions (9) and (10) (see A).

The solution for the Gaussian model (given by (8) with $\alpha=1$) is obtained from equation (11) taking $\alpha=1$,

[TABLE]

A very commonly used expression in the literature for the Gaussian model is obtained from the solution of the advection-diffusion equation (Eq. (8) with $\alpha=1$) in an infinite medium. For this equation to satisfy the boundary conditions given by equations (9) and (10), with $-\infty<x<\infty$ and $-\infty<z<\infty$, a ’mirror-image’ source is considered. The solution for this so called operational Gaussian model (O-G model) is [23]

[TABLE]

In the next section we are going to compare the solution (11) of our fractional model (8) against experimental data and both the Gaussian (13) and O-G (14) solution models.

3 Results and discussion

In order to analyze the performance of the fractional model proposed in this work, in comparison to the Gaussian model derived from Fick’s law, we consider only the simplest case of constant values for wind speed and eddy diffusivity. The average wind speed is obtained directly from the experimental data. To obtain a constant eddy diffusivity $K_{z}$, we follow the procedure introduced in [17], and we consider a spatial average $K_{z}=\langle K\rangle$ of a eddy diffusivity that is a linear function of downwind distance expressed by $K=\rho ux$, where $\rho$ is the turbulence parameter. The turbulence parameter $\rho$ is parameterized as the square of turbulent intensity using Taylor statistical theory of diffusion $\rho=(\frac{\sigma_{w}}{u})^{2}$ [26], where $\sigma_{w}$ is the standard deviation of the vertical wind speed component. In addition, to investigate the potential of application of fractional differential equations to model real data, we separated the experimental data from the Copenhagen experiment [18], the Prairie Grass experiment [19] and the Hanford experiment [20] into two groups: one with $\frac{h}{|L|}<10$ and another with $\frac{h}{|L|}>10$, where $L$ is the Monin-Obukhov length [27] and $h$ is height of PBL. The parameter $\frac{h}{|L|}$ is obtained from the energy balance equation in a turbulent flow. This parameter can be used to evaluate some characteristics of the physical structure of the turbulent flow. For $\frac{h}{|L|}<10$ we have a predominance of mechanical energy input (wind shear) in the turbulent flow. For $\frac{h}{|L|}>10$ we have a predominance of energy input by thermal convection in the turbulent flow. The analysis of these two situations aims to show that the value of $\alpha$ (order of the fractional equation) that best describes the concentration distribution, when confronting the proposed fractional model with the experimental data, is related to the physical structure of the turbulent flow. Usually, the performance of dispersion models is evaluated from a well know set statistical indices described by Hanna [28] 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.

3.1 Results for $\frac{h}{|L|}<10$

In order to estimate the better $\alpha$ value for each experiment, we analyzed the solution of our model from $\alpha=0.60$ to $\alpha=0.99$ by steps of $0.01$. An important result we found is that, regardless of the experiment, for $\frac{h}{|L|}<10$ the model (with constant velocity and constant eddy diffusivity) describes relatively well all experiments with $\alpha=0.72$.

In this case, for the Copenhagen experiment, Table 1 shows that, when compared to the Gaussian and Operational Gaussian (O-G) models using an average wind speed and eddy diffusivity, our model generates good results for the concentration distribution of contaminants in the PBL generated by a turbulent flow where the source of thermal convection and mechanical input is quite relevant. It also presents good results when compared to some Eulerian dispersion models found in the literature, which employ integer-order derivatives in the advection-diffusion equation (traditional integer-order models are derived from Fick’s law). In the case of the Moreira model [29], the stationary advection-diffusion equation employs a wind speed that is a function of height $z$, and a eddy diffusivity that is a function of height $z$ and horizontal distance $x$. The differential equation is solved by the GILTT method [30], extended to the case where the eddy diffusivity is a function of $z$ and $x$. In the case of the Kumar model [31] the GILTT method is used to solve the advection-diffusion equation, but in this case, the wind speed and eddy diffusivity are functions only of the height $z$. Table 1 shows that the fractional model proposed in this work, where an average wind speed and eddy diffusivity (therefore constant) was used, has a similar performance to the model of [29]. We also see that the fractional model proposed in this work performs better than Kumar model [31], which is also a traditional model (derived from the integer-order), but considers a wind speed and a eddy diffusivity dependent on the height $z$. We can also observe that a wind speed and eddy diffusivity that more correctly describes the turbulent flow, used in the models of Moreira [29] and Kumar [31], tends to compensate the deficiency of the mathematical structure of the classical advection-diffusion equation to describe the concentration distribution. In our work, a constant wind speed and eddy diffusivity were used precisely to show the ability of the fractional advection-diffusion equation to more accurately describe the concentration distribution of contaminants. In addition to the statistical indices, Figure 1 shows the scatter diagram of observed concentration and predicted concentration. Lines indicate a factor of two.

In Table 2 we show the results obtained for the stable Prairie Grass experiment with $\frac{h}{|L|}<4$. In this case, the stability regime can be considered closer to neutral. We can observe that the fractional model with $\alpha=0.72$ performs better than both Gaussian and O-G models for all distances. In addition, we see that for greater distances the difference between the fractional model and the Gaussian model is more evident. The Figure 2 shows the scatter diagram of observed concentration and predicted concentration for 200m and 800m distances. We can observe that the performance of the fractional model is superior to the performance of the Gaussian model at both distances. We attribute the better performance of the fractional model in relation to the Gaussian model to the fact that the dispersion of contaminants in a turbulent flow does not obey a Gaussian probability distribution but rather an anomalous diffusion distribution, induced by the turbulence, with power law mean square displacement [1].

Table 3 shows the results of the fractional, Gaussian and O-G models for the Hanford experiment with $\frac{h}{|L|}<10$. Again, the fractional model has a better performance than both Gaussian and O-G models. However, the difference between the fractional and the Gaussian model is smaller than in the cases of the Copenhagen and Prairie Grass experiments. The experiment of Copenhagen took place in an unstable regime where the source of thermal convection and mechanical input is quite relevant. In this case, the flow is strongly turbulent. The Prairie Grass experiment occurred in a stable regime, but with a high wind speed ($\approx 5.3ms^{-1}$). In this case, the mechanical source of the turbulence is strong. Furthermore, in the Prairie Grass experiment, small distances from the source are considered. It is reasonable to suppose that a diffusion model with constant wind speed and constant eddy diffusivity, that is, a model that does not completely describe the physical structure of the turbulent flow, performs better in a turbulent flow where the properties vary little with vertical height and for large distances from the source, where the dispersion of the contaminant can be described by the statistical properties of the turbulent flow. In the Hanford experiment, we have a stable regime, with a very low wind speed ($\approx 1.5ms^{-1}$). Therefore, the flow is slightly turbulent and the characteristics of anomalous diffusion, present in a strongly turbulent flow, are not very evident. In this case, it is expected that the Gaussian probability distribution can be employed to model the problem. Figure 3 shows the scatter diagram of observed concentration and predicted concentration for the Hanford experiment. We can see that the fractional and Gaussian models perform very similarly. However, the O-G model, given by equation (14), fails completely in this case. This is related to the approximation made so that the solution satisfies the boundary conditions (9) and (10). The values of the distances in the $z$ and $x$ directions involved in the experiments of Copenhagen and Prairie Grass are similar in almost all experiments. In this case, the approximation of boundary at infinity gives a good result. In particular, for the stable Prairie Grass where the $z$ and $x$ lengths are $\approx 400m$ and $800m$, respectively, the difference between Gaussian and O-G is negligible. In the Hanford experiment, the values of the distances in the $z$ and $x$ directions are very different ($\approx 200m$ and $3200m$, respectively), making the boundary at infinity approach totally inadequate.

3.2 Results for $\frac{h}{|L|}>10$

For $\frac{h}{|L|}>10$, Tables 4 and 5 show that the proposed fractional model with $\alpha=0.80$ performs better when compared to the Gaussian model in most statistical indices. As in the previous case, we analyzed the solutions of our model from $\alpha=0.60$ to $\alpha=0.99$ by steps of $0.01$, and we found that the model (with constant wind speed and constant eddy diffusivity) describes relatively well all experiments with $\alpha=0.80$ when $\frac{h}{|L|}>10$. We also observed that the performance of all models with $\frac{h}{|L|}>10$ is worst than the performance when $\frac{h}{|L|}<10$ (Tables 1, 2 and 3). This may be related to the difference in the structure of the turbulent flow in the PBL in both cases. When $\frac{h}{|L|}>10$ a free-convection-like state emerges [32]. In this case, we have a strong turbulence generated by thermal convection and a great variation in the structure of the flow in the vertical direction. This results in a large variation with the height of the intensity of the vertical eddy diffusivity. As in the experiments of Prairie Grass (unstable), and cases 1, 3, 7 and 8 of the Copenhagen experiment. On the other hand, when $\frac{h}{|L|}$ is very small ($\frac{h}{|L|}<10$) there is a turbulence where the mechanical source (wind shear) is relevant. In this case, we have a greater spatial homogenization in the flow. Consequently, we have a low variation with the height of the intensity of the vertical eddy diffusivity [33], as in the Prairie Grass (stable), and the Hanford experiments, and also in cases 2, 4, 5, 6 and 9 of the Copenhagen experiment. Models that employ a constant wind speed and eddy diffusivity will have greater difficulty in correctly describing the distribution concentration when there is a greater spatial asymmetry, as in the case of the PBL flow if $\frac{h}{|L|}>10$. This deficiency will be naturally compensated when non-constant wind speed and eddy diffusivity are used to more accurately describe the physical structure of the flow. However, in the present work we consider only constant wind speed and eddy diffusivity since, at this moment, we are interested only in demonstrating that the structure of the fractional order differential equation is more appropriate to calculate the distribution of dispersed contaminants in a turbulent flow than an integer-order differential equation. We are confronting a Gaussian model with fractional order derivatives against the Gaussian model with integer-order derivatives. Of course, a fractional or an integer-order model which considers a non-constant wind speed and eddy diffusivity, that more adequately describes the structure of the turbulent flow, will generate more realistic concentration distribution values. However, in general, it is not simple to obtain an analytical solution for both equations (fractional and integer-order) to easily compare its performances.

Finally, the most important result we found by considering both $\frac{h}{|L|}<10$ and $\frac{h}{|L|}>10$ is that there should be a relation between the order of the fractional derivative with the physical structure of the turbulent flow since, regardless of the experiment, when we have a predominance of mechanical energy input ($\frac{h}{|L|}<10$) in the turbulent flow all experiments where best described with $\alpha=0.72$, and when we have a predominance of energy input by thermal convection ($\frac{h}{|L|}>10$) the experimental data is better described by $\alpha=0.80$.

4 Conclusions

In the present work, we propose a new parameterization for the flow of concentration using fractional derivatives. The fractional order differential equation in the longitudinal and vertical directions is used to obtain the concentration distribution of contaminants in the Planetary Boundary Layer. The use of fractional derivatives in modeling diffusion of contaminants is justified by the presence of anomalous diffusion generated from 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 dealt 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 context, the use of fractional differential operators emerged as a valuable mathematical tool to model the time evolution of anomalous diffusion [1, 2, 4]. However, the use of fractional calculus to study steady-state regimes of dispersion of contaminants in the atmosphere is actually underexplored. In this respect, we recently demonstrated in [17] that the steady-state regime for a spatial concentration distribution of a non-reactive contaminants, in an anomalous diffusion process displaying a power law mean square displacement, is given naturally by a differential equation containing fractional derivatives in the advective term ($x$ direction) of the advection-diffusion equation. In the present work we propose a new parameterization for the flux of concentration using fractional derivatives (Eq. (5)). In this case the resulting fractional advection-diffusion equation containing fractional derivatives in the advective and diffusive terms ($x$ and $z$ directions), with this the model consider anomalous diffusion in both directions, differently of our previous fractional model [17] where the anomalous diffusion was considered only in the longitudinal direction. We solve the model and we compare the solution with both real experiments and traditional integer order derivative models. We show that our fractional model, even with constant wind speed and a constant eddy diffusivity, gives very good result in fitting the experimental data. The results obtained show that the structure of the fractional order differential equation is more appropriate to calculate the distribution of dispersed contaminants in a turbulent flow than an integer-order differential equation. Furthermore, a very important result we found it is that there should be a relation between the order $\alpha$ of the fractional derivative with the physical structure of the turbulent flow since, regardless the experiment, when we have a predominance of mechanical energy input in the turbulent flow all experiments where best described with $\alpha=0.72$, and when we have a predominance of energy input by thermal convection the experimental data is better described by $\alpha=0.80$. We are currently working towards a mathematical relationship between some parameters that describe the physical structure of the flow with the value of $\alpha$. Finally, the result obtained from the fractional derivative model motivates further investigations of applications of fractional differential equations to model diffusion of contaminants. There are many directions of investigation left to explore. Interesting examples include a theoretical investigation of a direct relationship between the fractal structure of turbulent flow and the power law mean square displacement for the anomalous diffusion (namely, a direct relation between the fractal structure and the order $\alpha$ of the fractional derivative), the use of local fractional/fractal derivatives in the model instead of the non-local Caputo derivatives, a detailed analyzes of the spatial memory effect introduced by a non-local derivative in the model, etc. These and other examples are left to future works.

Acknowledgments

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

Appendix A Solution of the fractional model

In this appendix we obtain the solution (11) of the fractional differential equation (8) subjected to the boundary conditions (9) and (10). The fractional differential equation (8) can be analytically solved by using the separation of variables technique. We begin by assuming that the solution can be written as follows,

[TABLE]

By substituting (15) into (8), we get two ordinary fractional differential equations in the variables $X$ and $Z$ as follows,

[TABLE]

and

[TABLE]

The solution for the system of equations (16) and (17) can be obtained by Frobenius method. The Eq. (16) has solution $X(x)=E_{\alpha}(-\kappa\lambda^{2}x^{\alpha})$, [17]. In order to solve (17) by Frobenius method, we consider the following power series

[TABLE]

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

[TABLE]

By considering $\beta=p+1$ we get

[TABLE]

where (Goulart et al., 2017),

[TABLE]

Then we have

[TABLE]

where, we use the identity $\Gamma(x+1)=x\Gamma(x)$.

From (21) we have two cases when $0<p\leq 1$ and $p=0$, respectively.

Case $0<p\leq 1$.

For the case $0<p\leq 1$ we have from (22)

[TABLE]

By substituting equation (23) into (17) we obtain

[TABLE]

that yields $p=\alpha$ and $b_{j+1}+\lambda^{2}b_{j}=0$. Then the solution of (17) when we have $0<p\leq 1$ is given by

[TABLE]

But this solution does not satisfy the boundary condition ${{}_{0}^{C}D^{\alpha}_{z}}Z_{1}(0)=0$ because

[TABLE]

then, if $b_{0}\neq 0$,

[TABLE]

Therefore $Z_{1}(z)$ does not gives the solution for our problem.

Case $p=0$.

For the case $p=0$ we have from (22),

[TABLE]

By substituting equation (28) into (17) we obtain

[TABLE]

that yields $b_{j+1}+\lambda^{2}b_{j}=0$. Then, the solution of (17) in this case is given by

[TABLE]

that gives the solution for our problem. From (30) we have

[TABLE]

By applying the boundary conditions (23) we get

[TABLE]

and

[TABLE]

Consequently, whe should have $\lambda=\lambda_{n}$, where $\lambda_{n}$ is a solution for the equation

[TABLE]

Finally, from the superposition principle we obtain (11):

[TABLE]

In order to fix the constants $a_{n}$ we use the boundary condition $c^{y}(0,z)=\frac{Q}{u}\delta(z-H_{s})$ in (11). We consider

[TABLE]

as a good approximation because, for $z\neq 0$, $E_{\alpha+1}(-\lambda_{n}^{2}z^{\alpha+1})$ goes to zero when $n$ goes to infinite. Now, by making

[TABLE]

where $p\in\{0,1,2,...,m\}$, we fix $a_{n}$ from the numerical solution of the following system:

[TABLE]