Simulation of atmospheric pollutant dispersion considering a bi-flux process and fractional derivatives

Atmospheric Pollution Research xxx (xxxx) xxx–xxx HOSTED BY

Contents lists available at ScienceDirect

Atmospheric Pollution Research journal homepage: www.elsevier.com/locate/apr

Original Article

Simulation of atmospheric pollutant dispersion considering a bi-flux process and fractional derivatives Anderson Palmeira, Paulo Xavier, Davidson Moreira∗ Manufacturing and Technology Integrated Campus, SENAI CIMATEC, Salvador, BA, Brazil

A R T I C LE I N FO

A B S T R A C T

Keywords: Laplace decomposition method Advection–diffusion equation Fractional derivative Bi-flow evolutionary model Anomalous dispersion

In this paper, we present an analytical solution of the fractional two-dimensional advection–diffusion equation with a bi-flow evolutionary model, which represents a modification of Fick's law applied to the dispersion of pollutants in the planetary boundary layer. The solution is obtained using Laplace decomposition to derive the Mittag–Leffler function, which is intrinsic to the solution of fractional differential equations. The solutions exhibit fast convergence and are in good agreement with data from the traditional Copenhagen (moderately unstable) and Hanford (stable to neutral) experiments in terms of the influence of the fractional derivative and the fourth-order term representing the retention phenomenon in the bi-flow evolutionary model. For the Copenhagen experiment, the best results are achieved with parameter values of α = 0.95 (fractional order of derivative) and β = 0.50 (retention control), with a normalized mean square error of 0.10 and 100% of results within a factor of two of the observed values. For the Hanford experiment, values of α = 1.00 and β = 1.00 are optimal, giving a normalized mean square error of 0.15 and 83% of results within a factor of two of the observed values, and exhibiting some dependence on the atmospheric stability.

1. Introduction The processes governing the transport and diffusion of pollutants in the atmosphere are numerous and very complex, and it is impracticable to describe them without the aid of mathematical models, especially the advection–diffusion equation. Such models allow the description, interpretation, and management of accidental releases, the evaluation of risk areas, identification of polluting sources, and the correction and control of pollutant concentrations in the atmosphere (Moreira and Vilhena, 2009). Although the advection–diffusion equation can be linear, it is still a challenging problem, because there are no analytical solutions when the physical parameters involved are dependent on all spatial and time variables, similar to the traditional turbulence closure problem. In this sense, the most common approach to solving the closure problem is based on the hypothesis of gradient transport (K theory or Fick's law). Analogously to molecular diffusion, this technique assumes that the turbulent flow of concentration is proportional to the magnitude of the mean concentration gradient (Mangia et al., 2002). Another important point is that the solutions found in the literature are usually for integer-order differential equations. Recently, there has been a resurgence in the study of fractional derivatives (Ross, 1974), which are a generalization of differential equations (Debnath, 2003;

Gorenflo and Mainardi, 2009; Moreira and Moret, 2018). The dispersion of pollutants is subject to atmospheric turbulence and, as a consequence, the phenomenon of anomalous diffusion arises naturally. Unlike ordinary diffusion, where the mean square displacement increases linearly over time, anomalous diffusion does not produce a linear mean square displacement. To deal with this phenomenon, two methods are widely used: the introduction of fractional derivatives to the model or the modification of Fick's law (Moreira et al., 2005a; Vilhena et al., 2008). Traditional differential equations (integer order) do not adequately describe the problem of turbulent diffusion, as the usual derivatives are not well defined in the non-differentiable regions introduced by turbulence. Consequently, the classical advection–diffusion equations are not expected to fully explain the anomalous diffusion of pollutants, because the system parameters generally grow faster than the solutions obtained by the classical models. In this way, fractional differential equations represent extensions of the basic equations of mathematical physics; they were explicitly introduced into the study of physical systems by Nigmatullin (1986) to describe anomalous diffusion in fractal geometry media (special types of porous media). Concerning the anomalous diffusion related to the problem of turbulence closure, an alternative is the addition to the advection–diffusion equation of a recent modification to

Peer review under responsibility of Turkish National Committee for Air Pollution Research and Control. ∗ Corresponding author. E-mail address: [email protected] (D. Moreira). https://doi.org/10.1016/j.apr.2019.09.015 Received 25 June 2019; Received in revised form 24 September 2019; Accepted 24 September 2019 1309-1042/ © 2019 Turkish National Committee for Air Pollution Research and Control. Production and hosting by Elsevier B.V. All rights reserved.

Please cite this article as: Anderson Palmeira, Paulo Xavier and Davidson Moreira, Atmospheric Pollution Research, https://doi.org/10.1016/j.apr.2019.09.015

Atmospheric Pollution Research xxx (xxxx) xxx–xxx

A. Palmeira, et al.

L [Lu (x , t )] + L [Ru (x , t )] + L [Nu (x , t )] + L [h (x , t )]

Fick's law, which is called the bi-flow diffusion process (Bevilacqua et al., 2011). The bi-flow process consists of the simultaneous flow of two sets of scattered particles with two distinct velocities. The first set obeys Fick's classical law, the so-called primary flow, whereas the second, called the secondary flow, follows a new ‘retention’ law. This gives rise to a fourth-order equation that provides a consistent approach for considering the rotational energy, which would otherwise not be possible with classical theory (Bevilacqua et al., 2016). This idea arises from evidence that particle retention may occur in diffusion processes for some dispersing substances immersed in a particular supporting media (Atsumi, 2002; D'Angelo et al., 2003; Deleersnijder et al., 2006). The classical diffusion equation does not represent the retention phenomenon satisfactorily, as there is a new phenomenon involved in the process that cannot be characterized by simply manipulating the diffusivity parameters. The mathematical description of the problem of atmospheric dispersion is a subject of great interest in the scientific community. The analytical solutions found in the literature are basically for linear problems that obey the traditional Fick's law, and differential equations of whole order; representative studies include those of Rounds (1955), Demuth (1978), Nieuwstadt and Haan (1981), Tirabassi et al. (1987), Tirabassi (1989), Tirabassi and Rizza (1994), Moreira et al., 2005b, Moreira et al., 2005c, 2009, 2014), Sharan and Modani (2006), Essa et al. (2007), Perez-Guerrero et al. (2012), and Pimentel et al. (2014). However, the problem of atmospheric dispersion has only recently been approached using fractional derivatives (Goulart et al., 2017; Moreira and Moret, 2018; Xavier et al., 2019; Acioli et al., 2019; Moreira and Santos, 2019; Moreira et al., 2019). Focusing our attention in this direction, this paper reports an analytical series solution for the fractional two-dimensional (2D) advection–diffusion equation considering nonFickian closure to simulate pollutant dispersion in the atmosphere. The novelty of this work relies on the analytical solution of this equation using the Laplace decomposition method (LDM) (Adomian, 1994; Jafari and Gejji, 2006) to consider the anomalous dispersion in the fractional derivative and bi-flux term in the vertical direction, which is more important in atmospheric problems due to the influence of the terrestrial surface. The LDM approach is based on a Laplace transform technique employed within a numerical algorithm to solve nonlinear ordinary and partial differential equations. The method is very well suited to physical problems as it does not require unnecessary linearization, perturbation, and other restrictive methods and assumptions which may change the problem being solved, sometimes seriously. To reach this goal, the remainder of this paper is organized as follows. In Section 2, we report the derivation of the solution for the fractional steady 2D advection–diffusion equation; in Section 3, the numerical results given by the method are presented and compared with experimental data; finally, in Section 4 we present the conclusions to this study.

(2)

resulting in

L [u (x , t )] =

f (x ) g (x ) 1 1 + 2 + 2 L [h (x , t )] − 2 L [Ru (x , t )] s s s s 1 − 2 L [Nu (x , t )] s

(3)

The second step is to represent the solution through an infinite series: ∞

u (x , t ) =

∑ un (x , t )

(4)

n=0

The nonlinear operator can be written as: ∞

Nu (x , t ) =

∑ An ,

(5)

n=0

where An are the Adomian polynomials, which can be calculated as: ∞

An =

1 dn ⎡ ⎛ ⎞⎤ N ∑ λiui⎟ ⎥ , n = 0,1,2, ... n! dλn ⎢ ⎜ i = 0 ⎠ ⎦λ = 0 ⎣ ⎝

(6)

Substituting Eqs. (4) and (5) into Eq. (3), we have: ∞

f (x ) g (x ) 1 1 ⎡ ⎤ + 2 + 2 L [h (x , t )] − 2 L L ⎢ ∑ un (x , t ) ⎥ = s s s s ⎣ n=0 ⎦ ∞

∞

1 ⎤ ⎡ ⎤ ⎡ ⎢R ∑ un (x , t ) ⎥ − s 2 L ⎢ ∑ An ⎥ ⎣ n=0 ⎦ ⎦ ⎣ n=0

(7)

For simplicity, the terms derived from the initial conditions and the source term can be separated as follows:

K (x , s ) =

f (x ) g (x ) 1 + 2 + 2 L [h (x , t )] s s s

(8)

so that, by applying the inverse Laplace transform to Eq. (8), we obtain:

1 L −1 [K (x , s )] = K 0 (x , t ) = f (x ) + tg (x ) + L −1 ⎡ 2 L [h (x , t )] ⎤ ⎦ ⎣s

(9)

Thus,

u 0 (x , t ) = K 0 (x , t )

(10)

Recursively, the other terms of the series can be written as follows:

1 1 un + 1 (x , t ) = −L −1 ⎡ 2 L [Run (x , t )] + 2 L [An ] ⎤ s ⎣s ⎦

,

n≥0 (11)

For more details about this method, see the work of Adomian (1994). Obviously, the solution given by this methodology depends too much on the choice of K 0 (x , t ) .

2. Methods 2.2. Fractional modelling: laplace transform 2.1. Brief description of LDM At this point, it is important to mention that fractional calculus is a new tool which has only been used in recent years to model complex systems with long-term memory. Despite an intricate mathematical background, fractional calculus is a generalization of differentiation and integration to arbitrary non-integer orders. In this context, fractional derivatives are a new tool for modelling complex systems, especially in atmospheric problems. In this sense, various tools for solving fractional equations have been reported in the literature, where the well-known Laplace transform is frequently applied. The Laplace transform of the Caputo fractional derivative is a generalization of the Laplace transform of integer-order derivatives, and is very important because it allows the use of initial values of classical integer-order derivatives with known physical interpretations. For a better understanding of the historical aspects and difficulties of fractional calculus,

LDM reduces the computational load and improves the accuracy of the results. As an example, consider the following nonlinear, nonhomogeneous, second-order partial differential equation of integer order:

Lu (x , t ) + Ru (x , t ) + Nu (x , t ) = h (x , t )

(1)

∂2 , ∂t 2

where L = R is the linear operator, N is the nonlinear operator, and h (x , t ) is the source term. We have the following initial conditions:

u (x , 0) = f (x )

,

∂u (x , 0) = g (x ) ∂t

(1a)

Using LDM, the first step is to apply the Laplace transform L in the variable t (t → s ) to both sides of Eq. (1): 2

Atmospheric Pollution Research xxx (xxxx) xxx–xxx

A. Palmeira, et al.

parameter of fractional order, c is the crosswind-integrated concentration, and β controls the balance between retention and diffusion. For simplicity, we take the vertical eddy diffusivity K z1 as a constant and, in a similar manner, K z2 is a constant related to the retention term. For dimensional consistency and following Gomez-Aguilar et al. (2016), an auxiliary parameter ϕ is introduced as follows (Moreira and Santos, 2019):

refer to Podlubny (1999). The memory and hereditary effects of conventional fractional derivatives have found numerous applications in science and engineering (Ionescu et al., 2017). However, the Laplace transform of a fractional derivative has not been explored in terms of atmospheric problems. The Laplace transform of a fractional derivative given by Caputo's formula is:

d 1 dα → 1−α α dx ϕ dx

n−1

L {Dtα f (t )} = s αF (s ) −

∑ s α−k −1f (k) (0) n − 1 < α ≤ n k=0

(12)

This is true if the parameter ϕ has dimensions of length. Therefore, the fractional representation of Eq. (17) with the inclusion of the parameter ϕ is as follows:

where s is the transformed variable, F (s ) = L {f (t ); s} , Dtα f (t ) is the derivative of arbitrary real order α of a function f (t ) . Under natural conditions on this function, for α → n (where n is an integer), the Caputo derivative becomes the conventional nth derivative of the function f (t ) . For more details, see the work of Podlubny (1999). Therefore, we can write:

L {Dtα f (t )} = s αF (s ) − s α − 1f (0), 0 < α ≤ 1

u

{

}

d f (t ) = sF (s ) − f (0) dt

(14)

2.3. Bi-flux model Consider a 1D diffusion process consisting of N particles moving in a homogeneous, isotropic supporting medium according to Fick's law. If a fraction f (N ) = (1 − β ) of the diffusing particles are temporarily trapped along their trajectories due to some mechanical, physical, chemical, or physicochemical interaction with the medium, the governing equation must include a fourth-order differential term of the form (Bevilacqua et al., 2011):

∂ 4c ∂z 4

Kz

∂c =0 ∂z

,

z=h

Kz

∂c =0 ∂z

,

z=0

(19)

(19a)

(19b)

In addition, an initial condition with constant emission Q is assumed at the source:

uc (0, z ) = Qδ (z − Hs )

,

x=0

(19c)

where δ (z − Hs ) is the Dirac delta function and Hs is the source height. The Dirac delta function can be approximated by: ∞

(15)

δ (z − Hs ) =

Then,

∂ 2c

0 < α ≤ 1,

where u = 1 − α . This fraction affects the advective term, increasing or ϕ decreasing the intensity of the wind speed. In typical physical situations with low levels of fractality, i.e. the order α of the fractional derivatives of the corresponding terms of the dynamic equations deviates slightly from an integer value n (in our case, 1), the term ϕ1 − α does not have a significant influence on the wind speed u. Therefore, it is assumed that ϕ = 1 [L], making the left side of Eq. (19) dimensionally consistent. Thus, the proposed procedure yields the solution of the problem without loss of generality. To obtain the solution of Eq. (19), it is necessary to specify the boundary conditions. As usual, the conditions of zero flux apply on the top of the PBL and on the ground:

(13)

There are several definitions of fractional-order derivatives, such as those of Riemann–Liouville, Caputo, Riesz, Weyl, Grunwald–Letnikov, and others. Here, Caputo's definition is used because it incorporates the initial conditions of the function, which makes it appropriate for solving applied problems such as air pollution modelling.

− β (1 − β ) K

∂αc ∂ 2c ∂ 4c = βK z1 2 − β (1 − β ) K z 2 4 , α ∂x ∂z ∂z U

If α = 1, Eq. (13) reduces to the classical Laplace transform of an integer-order derivative:

L

(18)

1⎡ ⎤ 1 + 2 ∑ cos(λn z )cos(λn Hs ) ⎥ h⎢ n=1 ⎦ ⎣

(20)

nπ h

∂ 4c

∂c = βK 2 − β (1 − β ) R 4 ∂t ∂z ∂z

are the eigenvalues. Therefore, the source condition (19c) and λn = can be rewritten as:

(16)

∞

The fourth-order term with the negative sign introduces the retention effect. The coefficients K and R are generalized constants. In atmospheric problems, K is the traditional eddy diffusivity [L2/T], whereas R is associated with the retention term [L4/T], which must be less than K. It is important to keep the parameters β and β (1 − β ) explicitly in the equation, because they express the balance between diffusion and retention when both are activated simultaneously. For more details about the proposed method, see Bevilacqua et al. (2011, 2016).

c (0, z ) =

Q ⎡ ⎤ 1 + 2 ∑ cos(λn z )cos(λn Hs ) ⎥ uh ⎢ n=1 ⎣ ⎦

(21)

For the solution of Eq. (19), the Laplace transform is applied to the variable x using the definition of the Caputo derivative, which gives:

∂ 2c ∂ 4c us αc‾ (s, z ) − us α − 1c (0, z ) = L ⎡βK z1 2 ⎤ − L ⎡β (1 − β ) K z 2 4 ⎤ ⎥ ⎢ ⎥ ⎢ ∂z ⎦ ∂z ⎦ ⎣ ⎣ (22) By rearranging Eq. (22), we obtain:

2.4. Fractional advection–diffusion equation with bi-flux term

c‾ (s, z ) = The fractional advection–diffusion equation including the bi-flow process can be written as follows:

∂αc ∂ 2c ∂ 4c U α = βK z1 2 − β (1 − β ) K z 2 4 , ∂x ∂z ∂z

0<α≤1

βK z1 ⎡ ∂2c ⎤ β (1 − β ) K z 2 ⎡ ∂ 4c ⎤ c (0, z ) + L L − 2⎥ 4 ⎢ ⎢ s us α ∂ z us α ⎦ ⎣ ⎦ ⎣ ∂z ⎥

(23)

Applying the inverse Laplace transform to Eq. (23) yields:

L −1 [c‾ (s, z )] = L −1 ⎡ ⎣

(17)

for 0 < z < h, x > 0 , where h is the planetary boundary layer height (PBL), x is the longitudinal source distance, z is the variable in the vertical direction, U is the longitudinal mean wind speed, α is the

βK z1 c (0, z ) ⎤ ∂ 2c + L −1 ⎡ α L ⎡ 2 ⎤ ⎤ ⎢ ⎢ s ⎦ ⎣ ∂z ⎥ ⎦⎥ ⎦ ⎣ us

− L −1 ⎡ ⎢ ⎣ 3

β (1 − β ) K z 2 ⎡ ∂ 4c ⎤ ⎤ L 4 ⎥ ⎢ us α ⎣ ∂z ⎥ ⎦⎦

(24)

Atmospheric Pollution Research xxx (xxxx) xxx–xxx

A. Palmeira, et al. 1

Then,

1

z 3 z 3 z z K z = 0.22w∗ h ⎛ ⎞ ⎛1 − ⎞ ⎡1 − exp ⎛−4 ⎞ − 0.0003 exp ⎛8 ⎞ ⎤ h⎠ ⎣ ⎝h⎠ ⎝ ⎝ h⎠ ⎝ h ⎠⎦

βK z1 β (1 − β ) K z 2 ⎡ ∂ 4c ⎤ ⎤ ∂ 2c c (x , z ) = c0 + L −1 ⎡ α L ⎡ 2 ⎤ ⎤ − L −1 ⎡ L ⎢ us ⎥ ⎢ ⎢ ⎥ ⎢ ∂z 4 ⎦ ⎥⎥ z us α ∂ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ (25)

(34) where w∗ is the convective velocity. For stable conditions, the vertical eddy diffusivity given by Ulke and Andrade (2001) can be used:

Following this methodology, the general solution can be written as: ∞

c (x , z ) = c0 + c1 + c2 + c3 + ...= ∑ ci i=0

(26)

z −1 K z = ku∗ z ⎛0.74 + 4.7 ⎞ L⎠ ⎝

(27)

In the model given by Eq. (33), the constant term K z1 can be obtained by averaging from the ground (z = 0) to the height h of the PBL. In fact, this average value is an approximation that limits the model, because it assumes homogeneous turbulence in the vertical direction (constant coefficients). However, the methodology proposed in this work allows a solution with variable coefficients. Here, for simplicity, the simulations use an average value for the vertical eddy diffusivities.

where ∞

c0 =

Q ⎡ ⎤ 1 + 2 ∑ cos(λn z )cos(λn Hs ) ⎥ uh ⎢ n=1 ⎣ ⎦

and

βK z1 β (1 − β ) K z 2 ⎡ ∂ 4ci ⎤ ⎤ ∂ 2c ci + 1 (x , z ) = L −1 ⎡ α L ⎡ 2i ⎤ ⎤ − L −1 ⎡ L ; 4 ⎥ ⎢ us ⎥ ⎢ ⎢ ⎥ ⎢ z us α ∂ ⎣ ⎦⎦ ⎣ ∂z ⎥ ⎦⎦ ⎣ ⎣ (28) i = 0,1,2,3, ...

3. Numerical results To evaluate the performance of the model derived in this study, simulation results were compared with the data observed in the traditional Copenhagen and Hanford experiments, which are used to obtain the crosswind-integrated concentrations.

For simplicity, we set A1 = βK z1 and A2 = β (1 − β ) K z 2 , resulting in:

c1 = −

c2 =

2Q xα uh uΓ (α + 1)

∞

∑ λn2 cos(λn z )cos(λn Hs )(A1 + A2 λn2) n=1

2Q x 2α uh 2!u2Γ (2α + 1)

c3 = −

(29)

∞

∑ λn4 cos(λn z )cos(λn Hs )(A1 + A2 λn2)2 n=1

2Q x 3α 3 uh 3!u Γ (3α + 1)

3.1. Copenhagen experiment (30) The Copenhagen dispersion experiments, described by Gryning and Lyck (1984) and Gryning et al. (1987) and then updated by Gryning and Lyck (2002), involved the release of SF6 tracer (sulphur hexafluoride) in the north of Copenhagen. The tracer was released without buoyancy from a tower with a height of 115 m and concentrations were measured at ground level. The sampling units were positioned at distances of 2–6 km from the release point. The surface roughness length was 0.6 m. The parameters describing the meteorological conditions and concentrations during the experiments are listed in Table 1. Note that the concentration was normalized by the emission rate (c / Q) . The parameters U10 and U115 indicate the wind speed measured at heights of 10 m and 115 m, respectively.

∞

∑ λn6 cos(λn z )cos(λn Hs )(A1 + A2 λn2)3 n=1

(31)

⋮ where Γ is the Gamma function. Therefore, grouping these first terms, we have:

c (x , z ) =

+

∑n = 1 cos(λn z )cos(λn Hs ) ⎧1 − ⎨ ⎩ 2 cos(λn z )cos(λn Hs )(A1 + A2 λn )+ Q uh

+

2Q uh

x 2α 1 λ4 2 ! u2Γ (2α + 1) n

∞

xα λ2 uΓ (α + 1) n

1

cos(λn z )cos(λn Hs )(A1 + A2 λn2)2 − 3 !

cos(λn z )cos(λn Hs )(A1 + A2 λn2)3 + ...⎫ ⎬ ⎭

x 3α λ6 u3Γ (3α + 1) n

Table 1 Micrometeorological data and crosswind-integrated concentrations observed during the Copenhagen diffusion experiments.

(32)

Thus, we obtain the following final solution:

c (x , z ) =

Q uh 2Q + uh

(35)

Run

U10 (ms−1)

U115 (ms−1)

h (m)

W

1

2.1

3.4

1980

1.8

2

4.9

10.6

1920

1.8

3

2.4

5.0

1120

1.3

4 5

2.5 3.1

4.6 6.7

390 820

0.7 0.7

6

7.2

13.2

1300

2.0

7

4.1

7.6

1850

2.2

8

4.2

9.4

810

2.2

9

5.1

10.5

2090

1.9

∞

∑

cos(λn z )cos(λn Hs )

n=1 xα

Eα ⎡− λn2 (βK z1 + β (1 − β ) K z 2 λn2) ⎤ ⎦ ⎣ u

(33)

Equation (33) is quite general in the sense that, if K z2 = 0 , β = 1, and α = 1, we have the solution for normal diffusion. Note that Eα is the Mittag–Leffler function, which is intrinsic to solutions with fractional derivatives. The solution given by Eq. (33) represents an important advance in analytical solutions of the fractional advection–diffusion equation using LDM. Although the parameters considered in this work are constants, the solution procedure can be used for coefficients that vary in the vertical direction. 2.5. Turbulence parameterisation In this work, for moderately unstable conditions, was assumed the vertical eddy diffusivity proposed by Degrazia et al. (1997): 4

⃗

(ms−1)

x (m)

c/Q (10−4 sm−2)

1900 3700 2100 4200 1900 3700 5400 4000 2100 4200 5100 2000 4200 5900 2000 4100 5300 1900 3600 5300 2100 4200 6000

6.28 2.31 5.38 2.95 8.2 6.22 4.3 11.66 6.72 5.84 4.97 3.96 2.22 1.83 6.7 3.25 2.23 4.16 2.02 1.52 4.58 3.11 2.59

Atmospheric Pollution Research xxx (xxxx) xxx–xxx

A. Palmeira, et al.

Fig. 1. Convergence test of the model with β = 0.50 and different values of α (0.96, 0.98, and 1.00).

As the proposed solution of the 2D advection–diffusion equation of fractional order is a series solution, a series convergence test is appropriate. For this test, we use meteorological data from Copenhagen experiment 1. Fig. 1 shows the numerical convergence of the proposed solution for the ground-level concentration at a distance of 1900 m as a function of the number of terms n (LDM solution) considering the wind speed U115 (measured at the height of the source), retention parameter β = 0.50 , and different values of the fractional parameter α (0.96, 0.98 and 1.00). Fig. 1 shows that the solution converges faster as the number of terms of the sum (n) increases, even for different values of the fractional parameter α . At a distance of 1900 m, a decrease in α increases the concentration. Fig. 2 shows the results of sensitivity tests to verify the influence of the fourth-order term, but only with an integer-order fractional parameter (α = 1.00 ). The vertical profile of the crosswind-integrated concentration at distances of 250 m and 2500 m from the source are presented for various values of β (0.10, 0.25, 0.50, 0.75, and 1.00). The simulations used an average value for Kz1 given by Eq. (34) and Kz2 = 1 m4/s. In Fig. 2(a), which represents a distance of 250 m, a decrease in the parameter β produces an increase in concentration in the region of the height of the source. Fig. 2(b), which represents a distance of 2500 m, exhibits a similar effect, but with a greater homogenization of the concentration in the vertical direction. The parameter β effectively controls the diffusive effect of pollutant dispersion in the atmosphere. Fig. 3 shows the crosswind-integrated concentration as a function of the source distance for various values of β and α . Fig. 3(a) shows that, by decreasing the value of β , the peak concentration can be shifted to greater distances from the source, but the peak concentration value does not change. However, in Fig. 3(b), we can see a slight decrease in the peak concentration when α decreases, but the concentration level further from the source increases. This is more pronounced at the lower value of β . It is important to demonstrate the influence of the parameter β on the ground-level concentration. This effect is shown in Fig. 4. It is clearly apparent that β effectively influences the crosswindintegrated concentration at the ground level. Therefore, the fraction β is a control parameter for the intensity of the secondary flux. The observed concentrations are compared with those from simulations with wind speed U115 U115 and α = 0.98 for various β in Fig. 5. The crosswind-integrated concentration with constant eddy

Fig. 2. Vertical profile of crosswind-integrated concentration for distances of a) 250 m and b) 2500 m from the source as a function of β (0.10, 0.25, 0.50, 0.75, and 1.00).

diffusivity modelled under the moderately unstable conditions of the Copenhagen experiment is very well represented using β = 0.50, with all of the data points lying between the dotted lines (representing a factor of two). The following analysis is a statistical evaluation of the model simulations in comparison to the experimental 2D data from the Copenhagen experiments. The aim is to verify the influence of the fractional parameter (α ) and retention term (β ), and to compare our results with those of other models. The model performance is statistically evaluated using the procedure described by Hanna (1989) and the following metrics: NMSE (normalized mean square error) = (Co − Cp)2 / Cp Co , FAT2 = fraction of data for which 0.5 ≤ (Cp/ Co) ≤ 2 , COR (correlation coefficient) = (Co − Co )(Cp − Cp / σo σp , FB (fractional bias) = Co − Cp /0.5(Co + Cp ) , FS (fractional standard deviation) = (σo − σp)/0.5(σo + σp),where the subscripts o and p refer to the observed and predicted quantities, respectively, and an overbar indicates an average value. The FB statistical index reflects whether the expected quantities underestimate or overestimate the observed values. The statistical NMSE index represents

5

Atmospheric Pollution Research xxx (xxxx) xxx–xxx

A. Palmeira, et al.

Fig. 4. Ground-level concentration as a function of β .

Fig. 3. Ground-level crosswind-integrated concentration as a function of source distance: a) α = 1.00 and β = 0.10, 0.25, 0.50, 0.75, and 1.00; b) α = 0.96, 0.98, and 1.00 and β = 0.5, 0.75.

the dispersion of the model output in relation to the dispersion of data. The best results are expected to give values close to zero for the NMSE, FB, and FS indexes, and close to 1 for the COR and FAT2 indexes. Note that the Copenhagen experiment provides wind speed data at two distinct heights: 10 m (U10 ) and 115 m (U115 ). Therefore, it is possible to analyse the solution obtained in this work and the influence of the height of wind speed measurements in the dispersion process. Table 2 presents the statistical results of the model represented by Eq. (33) for different values of α and β and for wind speed measured at heights of 10 m and 115 m. In addition, we present the results given by other models: Model 2 refers to the ADMM methodology (Advection–Diffusion Multilayer Technique) with eddy diffusivities dependent on x and z (inhomogeneous turbulence in the vertical direction), but using a mean eddy diffusivity in the x direction (Costa et al., 2006), and a wind speed U with a similar profile; Model 3 refers to the GILTT methodology (Generalized Integral Laplace Transform Technique) (Moreira et al., 2009), also with eddy diffusivity dependent on the vertical direction z, but with a wind speed profile U given by a power law. Regarding the model proposed in this work, Table 2 indicates that

Fig. 5. Observed and predicted scatter diagram of ground-level concentrations with α = 0.98 and U115 for different values of β (0.10, 0.25, 0.50, 0.75, and 1.00). Dotted lines indicate a factor of two.

the best result in comparison with the other models is achieved with parameter values of α = 0.98, β = 0.50, and wind speed U115, as this gives the lowest NMSE (0.10) and a FAT2 value of 100% (1.00). The result with α = 0.99 and β = 0.50 is very similar. The worst result occurs with α = 0.98 and wind speed U10, as this has the highest NMSE (0.81) and the lowest FAT2 (0.26). The fact that the physical parameters do not depend explicitly on the variable z, which is influenced by the earth's surface, does not affect the proposed methodology, as this lack of inhomogeneous turbulence in the vertical direction is somehow compensated by the variation of α and β . This is confirmed by the comparison between the proposed model and the results of Models 2 and 3 (ADMM and GILTT, respectively), which consider inhomogeneity in the vertical direction and are more complex models with respect to their development and computational implementation. The model proposed in this study is simpler and easier to implement, and allows 6

Atmospheric Pollution Research xxx (xxxx) xxx–xxx

A. Palmeira, et al.

Table 2 Statistical evaluation for Copenhagen experiment.

Model 1 (U115)

Model 1 (U10)

Model 2 (U similarity) Model 3 (U power law)

Model 1 (U115)

Model 1 (U10)

Model 1 (U115)

Model 1 (U10)

β = 0.10 β = 0.25 β = 0.50 β = 0.75 β = 1.00 β = 0.10 β = 0.25 β = 0.50 β = 0.75 β = 1.00 ADMM (Costa et al., 2006) GILTT (Moreira et al., 2009) β = 0.10 β = 0.25 β = 0.50 β = 0.75 β = 1.00 β = 0.10 β = 0.25 β = 0.50 β = 0.75 β = 1.00 β = 0.10 β = 0.25 β = 0.50 β = 0.75 β = 1.00 β = 0.10 β = 0.25 β = 0.50 β = 0.75 β = 1.00

NMSE

COR

FAT2

FB

FS

0.32 0.10 0.12 0.21 0.32 0.76 0.40 0.17 0.10 0.09 0.06

0.32 0.81 0.87 0.87 0.87 0.73 0.84 0.85 0.85 0.85 0.89

0.78 0.91 0.96 0.91 0.78 0.30 0.65 0.87 0.91 0.96 1.00

−0.11 −0.03 0.17 0.32 0.43 −0.73 −0.52 −0.26 −0.09 0.04 0.02

0.12 0.24 0.28 0.36 0.43 −0.41 −0.40 −0.25 −0.12 −0.03 0.09

0.09

0.85

1.00

0.11

0.13

0.36 0.10 0.11 0.18 0.28 0.79 0.43 0.19 0.11 0.09 0.40 0.12 0.10 0.17 0.25 0.81 0.46 0.22 0.12 0.09

0.26 0.79 0.86 0.87 0.87 0.69 0.84 0.85 0.85 0.85 0.20 0.76 0.86 0.87 0.87 0.65 0.84 0.85 0.85 0.85

0.78 0.87 1.00 0.91 0.83 0.30 0.65 0.89 0.91 0.96 0.74 0.87 1.00 0.91 0.87 0.26 0.61 0.87 0.87 0.91

0.11 0.04 0.14 0.28 0.39 −0.74 −0.54 −0.29 −0.13 0.00 0.09 0.05 0.12 0.25 0.36 −0.75 −0.56 −0.33 −0.16 −0.04

0.10 0.25 0.27 0.34 0.41 −0.40 −0.40 −0.27 −0.14 −0.05 0.06 0.25 0.27 0.33 0.40 −0.39 −0.40 −0.28 −0.17 −0.08

Table 3 Micrometeorological data and crosswind-integrated concentrations observed during the Hanford diffusion experiments. Run

U2 (ms−1)

L (m)

h (m)

u∗ (ms−1)

x (m)

c/Q (10−3 sm−2)

1

3.63

165

325

0.40

2

1.42

44

135

0.26

3

2.02

77

182

0.27

4

1.50

34

104

0.20

5

1.41

59

157

0.26

6

1.54

71

185

0.30

100 200 800 1600 3200 100 200 800 1600 3200 100 200 800 1600 3200 100 200 800 1600 3200 100 200 800 1600 3200 100 200 800 1600 3200

19.5 11.7 3.73 2.14 1.30 51.9 36.7 12.90 9.08 7.22 27.1 18.1 5.91 3.31 1.79 91.8 48.6 20.10 13.10 9.15 83.9 42.4 10.50 8.61 6.64 88.4 61.1 13.40 6.15 3.11

decrease in β leads to an increase in concentration in the region of the height of the source. A similar effect can be observed in Fig. 6(b), which represents a distance of 2000 m, although there is a greater homogenization of the concentration in the vertical direction. Again, the parameter β effectively controls the diffusive effect in the pollutant dispersion. Fig. 7 shows the crosswind-integrated concentration as a function of distance from the source for various values of β and α . Fig. 7(a) shows that a decrease in the parameter β leads to the concentration remaining larger further from the source. However, in Fig. 7(b), which shows the effect of varying α with β = 0.50, a decrease in α causes slightly higher concentrations with increasing distance from the source. In addition, a comparison between the predicted concentrations using a fractional parameter of integer order and the observed concentrations is presented in Fig. 8. Fig. 8 shows that the crosswind-integrated concentration with constant eddy diffusivity modelled under the stable conditions of the Hanford experiment is best represented by β = 1.00, with almost all of the data points lying between the dotted lines (representing a factor of two). Table 4 presents the statistical results of the model represented by Eq. (33) for different values of α and β . Normally, the Hanford experiment is used to study the effect of deposition on the soil. However, data are also provided for the non-depositing pollutant SF6. The simulations only considered distances of 800, 1600, and 3200 m, because these are the distances that have measured concentration data for the depositing pollutant ZnS. Thus, in the work of Moreira et al. (2010), the ADMM (derivative of integer-order) methodology was used to simulate SF6 for the same distances as for the depositing pollutant. Regarding the model proposed in this work, Table 4 indicates that the best result in comparison with the ADMM model was achieved with α = 1.00 and β = 1.00 (considering only distances of 800, 1600, and 3200 m), as this gave the lowest NMSE (0.15) and the highest FAT2

explicit analysis of its physical parameters because of its analytical characteristics. 3.2. Hanford experiment The Hanford dispersion experiments were conducted over May–June 1983 in a semi-arid region of south-eastern Washington on generally flat terrain. A detailed description of the experiment is provided by Doran et al. (1984) and Doran and Horst (1985). Data were obtained from six dual-tracer releases at 100, 200, 800, 1600, and 3200 m from the source during moderately stable to near-neutral conditions. A total of six test runs were conducted (Table 3). The data collected during the field tests were tabulated (as crosswind-integrated tracer concentration data) and reported by Doran et al. (1984). The terrain roughness was 3 cm. Two tracers, one depositing and one nondepositing, were released simultaneously from a height of 2 m (collected at 1.5 m). Zinc sulphide (ZnS) was used as the depositing tracer, and sulphur hexafluoride (SF6) was the non-depositing tracer. In this study, we consider only the non-depositing SF6 tracer. The micrometeorological data are listed in Table 3, where L is the Monin–Obukhov length, u∗ is the friction velocity, and U2 is the wind speed measured at 2 m. The observed crosswind-integrated concentration data at various distances were normalized according to the release rate Q. For the stable Hanford case, Fig. 6 shows the simulated vertical profile of the crosswind-integrated concentration at distances of 100 m and 2000 m from the source for various values of β and α = 1.00. The meteorological parameters of Hanford experiment 1 were used, with the average value of Kz1 given by Eq. (35) and Kz2 = 1 m4/s. Fig. 6(a), which represents a distance of 100 m, shows that a 7

Atmospheric Pollution Research xxx (xxxx) xxx–xxx

A. Palmeira, et al.

Fig. 6. Vertical profile of crosswind-integrated concentration for distances of a) 100 m and b) 2000 m from the source as a function of β = 0.10, 0.25,0.50, 0.75, and 1.00 for Hanford experiment 1 (α = 1.00).

Fig. 7. Ground-level crosswind-integrated concentration as a function of source distance: a) α = 1.00 and β = 0.10, 0.25, 0.50, 0.75, and 1.00; b) α = 0.96, 0.98, and 1.00 and β = 0.50.

value of 83% (0.83). The worst result was produced with α = 0.97 and β = 1.00 (although still better than the ADMM results), as this combination gave the highest NMSE (0.26) and the lowest FAT2 (0.72). Considering all distances, the simulations with α = 1.00 and β = 1.00 produced slightly better results, as FAT2 = 83%. However, the results are very similar to those with β = 0.75 and α = 0.97, 0.98, and 0.99. Therefore, for the Hanford stable experiment, there were no major differences in results. At this point, it should be noted that with weak winds (U < 2 ms−1), the longitudinal diffusive term should be considered, but this is not taken into account in the proposed model (Moreira et al., 2005c). In this sense, the parameters α and β cannot efficiently circumvent this physical situation. Even so, particularly when compared to a more complex model such as ADMM, the results can be considered quite satisfactory.

simultaneously and combined with the concept of fractional derivatives. We showed that the memory effect, which is normally only considered when the eddy diffusivity is dependent on the source distance, is vital for both turbulence closure and fractional derivatives. This is an interesting result, because dispersion problems in the atmosphere generally only take account of the memory effect for the eddy diffusivity. The best simulation results for the Copenhagen experiment, in comparison with the other models, was achieved with α = 0.98, β = 0.50, and wind speed U115. These parameters gave the lowest NMSE (0.10) and the maximum FAT2 of 100% (1.00). However, simulations with α = 0.99 and β = 0.50 gave very similar results. The worst result occurred with α = 0.98 and wind speed U10, giving the highest NMSE (0.81) and the lowest FAT2 (0.26). Regarding the Hanford experiment, the best result in comparison with the ADMM model was achieved using α = 1.00 and β = 1.00 (considering only distances of 800, 1600, and 3200 m), with a low NMSE (0.15) and FAT2 value of 83% (0.83). The worst result occurred with α = 0.97 and β = 1.00 (although this was still better than the results given by the ADMM model); in this case, the NMSE was 0.26 and the FAT2 was 0.72. Considering all distances, the

4. Conclusions This study has proposed the use of a new methodology to obtain solutions of the 2D advection–diffusion equation of fractional order to simulate the dispersion of pollutants emitted in the planetary boundary layer. To obtain this solution, LDM and bi-flux closure were applied 8

Atmospheric Pollution Research xxx (xxxx) xxx–xxx

A. Palmeira, et al.

parameterisations in the representation of physical processes. The fact that the physical parameters of the model (eddy diffusivity and wind speed) do not depend explicitly on the variable z, which is influenced by the earth's surface, does not affect the proposed methodology, as this lack of inhomogeneous turbulence in the vertical direction is somehow compensated by the variation of α and β . This is confirmed by the comparison between the model proposed in our study and the results of ADMM and GILTT, which consider the inhomogeneity in the vertical direction and are more complex models. The model proposed in this paper is simpler and easier to implement, and permits the explicit analysis of its physical parameters because of its analytical characteristics. In practical terms, we verified that there was no advantage in using an eddy diffusivity that is dependent on x and z simultaneously in the dispersion process, considering the meteorological conditions of the Copenhagen and Hanford experiments. Finally, this work is a step towards solving the advection–diffusion equation with fractional parameters in the diffusive terms and inhomogeneous turbulence in the vertical direction. This is a challenge, but our results show that new models should take into account the fundamental aspects of fractional modelling as a generalization of the calculations we currently know. Declaration of competing interest Fig. 8. Observed and predicted scatter diagram of ground-level concentrations with α = 1.00 for different values of β . Dotted lines indicate a factor of two.

We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Table 4 Statistical evaluation for Hanford experiment.

Model 1

800, 1600, 3200 m Model 2 800, 1600, 3200 m

Model 1

800, 1600, 3200 m

Model 1

800, 1600, 3200 m

Model 1

800, 1600, 3200 m

β = 0.10 β = 0.25 β = 0.50 β = 0.75 β = 1.00 β = 1.00 ADMM (Moreira et al., 2010) β β β β β β β β β β β β β β β β β β

= 0.10 = 0.25 = 0.50 = 0.75 = 1.00 = 1.00 = 0.10 = 0.25 = 0.50 = 0.75 = 1.00 = 1.00 = 0.10 = 0.25 = 0.50 = 0.75 = 1.00 = 1.00

Acknowledgement

NMSE

COR

FAT2

FB

FS

0.97 0.38 0.25 0.33 0.46 0.15 0.17

0.89 0.89 0.89 0.89 0.89 0.95 0.91

0.13 0.37 0.67 0.77 0.83 0.83 0.67

−0.79 −0.45 −0.15 0.04 0.17 −0.33 −0.32

−0.30 −0.02 0.25 0.41 0.53 −0.08 −0.10

1.02 0.41 0.25 0.31 0.43 0.18 1.07 0.44 0.25 0.29 0.40 0.22 1.13 0.48 0.26 0.28 0.37 0.26

0.89 0.89 0.89 0.89 0.89 0.95 0.89 0.89 0.89 0.89 0.89 0.95 0.89 0.89 0.89 0.89 0.89 0.95

0.13 0.33 0.67 0.80 0.80 0.78 0.13 0.33 0.60 0.80 0.77 0.72 0.13 0.30 0.50 0.80 0.77 0.72

−0.81 −0.48 −0.18 0.00 0.14 −0.37 −0.83 −0.51 −0.22 −0.03 0.10 −0.42 −0.85 −0.54 −0.25 −0.06 0.07 −0.46

−0.31 −0.04 0.23 0.39 0.51 −0.12 −0.32 −0.05 0.21 0.37 0.48 −0.16 −0.33 −0.07 0.19 0.35 0.46 −0.20

We are thankful for the support of SENAI CIMATEC and FAPESB, Brazil. References Acioli, P.S., Xavier, F.A., Moreira, D.M., 2019. Mathematical model using fractional derivatives applied to the dispersion of pollutants in the Planetary Boundary Layer. Boundary-Layer Meteorol. 170, 285–304. Adomian, G., 1994. Solving Frontier Problem of Physics: the Decomposition Method. Springer Netherlands, pp. 354. Atsumi, H., 2002. Hydrogen bulk retention in graphite and kinetics of diffusion. J. Nucl. Mater. 307 (311), 1466–1470. Bevilacqua, L., Galeão, A.C.N.R., Costa, F.P., 2011. A new analytical formulation of retention effects on particle diffusion processes. Annal. Brazilian Acad. Sci 83 (4), 1443–1464. Bevilacqua, L., Jiang, M., Silva Neto, A., Galeão, A.C.R.N., 2016. An evolutionary model of bi-flux diffusion processes. Braz. Soc. Mech. Sci. Eng. 38 (5), 1421–1432. Costa, C.P., Vilhena, M.T., Moreira, D.M., Tirabassi, T., 2006. Semi-analytical solution of the steady three-dimensional advection-diffusion equation in the PBL. Atmos. Environ. 40, 5659–5669. Debnath, L., 2003. Recent applications of fractional calculus to science and engineering. Int. J. Math. Math. Sci. 2003 (54), 3413–3442. Degrazia, G.A., Campos Velho, H.F., Carvalho, J.C., 1997. Nonlocal exchange coefficients for the convective boundary layer derived from spectral properties. Contrib. Atmos. Phys. 70, 57–64. Deleersnijder, E., Beckers, J.M., Delhez, E.M., 2006. The residence time of setting in the surface mixed layer. Environ. Fluid Mech. 6, 25–42. Demuth, C., 1978. A contribution to the analytical steady solution of the diffusion equation for line sources. Atmos. Environ. 12, 1255–1258. Doran, J.C., Horst, T.W., 1985. An evaluation of Gaussian plume-depletion models with dual-tracer field measurements. Atmos. Environ. 19, 939–951. Doran, J.C., Abbey, O.B., Buck, J.W., Glover, D.W., Horst, T.W., 1984. Field Validation of Exposure Assessment Models. Data Environmental Science Research Laboratory, U.S. Environmental Protection Agency, Research Triangle Park, NC, pp. 177pp EPA/600/ 384/092A. D'Angelo, M.V., Fontana, E., Chertcoff, R., Rosen, M., 2003. Retention phenomena in nonNewtonian fluids flow. Physics 327, 44–48. Essa, K.S.M., Etman, S.M., Embaby, M., 2007. New analytical solution of the dispersion equation. Atmos. Res. 84, 337–344. Gomez-Aguilar, J.F., Miranda-Hernandez, M., Lopez-Lopez, M.G., Alvarado-Martinez, V.M., Baleanu, D., 2016. Modeling and simulation of the fractional space–time diffusion equation. Commun. Nonlinear Sci. Numer. Simul. 30, 115–127. Gorenflo, R., Mainardi, F., 2009. Some recent advances in theory and simulation of fractional diffusion processes. J. Comput. Appl. Math. 229 (2), 400–415. Goulart, A.G.O., Lazo, M.J., Suarez, J.M.J., Moreira, D.M., 2017. Fractional derivative

simulations with α = 1.00 and β = 1.00 present slightly better results, as FAT2 = 83%. It should be noted that the proposed model with constant vertical eddy diffusivity and wind speed limit the problem of the dispersion of pollutants in the atmosphere. This approximation is responsible for discrepancies between the experimental data and the results of the proposed solution. However, in spite of this limitation, the methodology presented in this work can be used for variable coefficients in the vertical direction, without compromising the importance of fractional derivatives and its application with different turbulent 9

Atmospheric Pollution Research xxx (xxxx) xxx–xxx

A. Palmeira, et al.

6289–6294. Moreira, D.M., Xavier, P.H.F., Palmeira, A.S., Nascimento, E.S., 2019. New approach to solving the atmospheric pollutant dispersion equation using fractional derivatives. Int. J. Heat Mass Transf. 144, 118667. Nieuwstadt, F.T.M., Haan, B.J., 1981. An analytic solution of the one-dimensional diffusion equation in a non-stationary boundary layer with an application to the inversion rise fumigation. Atmos. Environ. 15, 845–851. Nigmatullin, R.R., 1986. The realization of the generalized transfer equation in a medium with fractal geometry. Phys. Status Solidi B 133, 425–430. Perez-Guerrero, J.S., Pimentel, L.C.G., Ulke, A.G., Oliveira-Junior, J.F., Heilbron Filho, P.F.L., 2012. A unified analytical solution of the steady-state atmospheric diffusion equation. Atmos. Environ. 55, 201–212. Pimentel, L.C.G., Perez-Grerrero, J.S., Ulke, A.G., Duda, F.P., Heilbron Filho, P.F.L., 2014. Assessment of the unified analytical solution of the steady-state atmospheric diffusion equation for stable conditions. Proc. R. Soc. A Math. Phys. Eng. Sci. 470 2014002120140021. Podlubny, I., 1999. Fractional Differential Equations. Academic Press, pp. 340. Ross, B., 1974. Fractional calculus and its applications. In: Proceedings of the International Conference. June. Springer Verlag, New Haven, pp. 386 New York. Rounds, W., 1955. Solutions of the two-dimensional diffusion equations. Trans. Am. Geophys. 36, 395–405. Sharan, M., Modani, M., 2006. A two-dimensional analytical model for the dispersion of air-pollutants in the atmosphere with a capping inversion. Atmos. Environ. 40, 3479–3489. Tirabassi, T., 1989. Analytical air pollution advection and diffusion models. Water Air Soil Pollut. 47, 19–24. Tirabassi, T., Rizza, U., 1994. Applied dispersion modelling for ground-level concentrations from elevated sources. Atmos. Environ. 28 (4), 611–615. Tirabassi, T., Tagliazucca, M., Galliani, G., 1987. Easy to use air pollution model for turbulent shear flow. Environ. Softw 2 (1), 37–44. Ulke, A.G., Andrade, M.F., 2001. Modeling urban air pollution in São Paulo, Brazil: sensitivity of model predicted concentrations to different turbulence parameterizations. Atmos. Environ. 35, 1747–1763. Vilhena, M.T., Costa, C.P., Moreira, D.M., Tirabassi, T., 2008. A semi-analytical solution of the three-dimensional advection-diffusion equation considering non-local turbulence closure. Atmos. Res. 90, 63–69. Xavier, P.H.F., Nascimento, E.G.S., Moreira, D.M., 2019. A model using fractional derivatives with vertical eddy diffusivity depending on the source distance applied to the dispersion of atmospheric pollutants. Pure Appl. Geophys. 176, 1797–1806.

models for atmospheric dispersion of pollutants. Physica A 477, 9–19. Gryning, S.E., Lyck, E., 1984. Atmospheric dispersion from elevated sources in an urban area: comparison between tracer experiments and model calculations. J. Clim. Appl. Meteorol. 23 (4), 651–660. Gryning, S.E., Lyck, E., 2002. The Copenhagen Tracer Experiments: Reporting of Measurements. Riso National Laboratory, pp. 74. Gryning, S.E., Holtslag, A.M.M., Irwin, J., Sivertsen, B., 1987. Applied dispersion modellingbased on meteorological scaling parameters. Atmos. Environ. 21, 79–89. Hanna, S.R., 1989. Confidence limit for air quality models as estimated by bootstrap and jacknife resampling methods. Atmos. Environ. 23, 1385–1395. Ionescu et al., 2017. Jafari, H., Gejji, V.D., 2006. Revised Adomian decomposition method for solving a system of nonlinear equations. Appl. Math. Comput. 175, 1–7. Mangia, C., Moreira, D.M., Schipa, I., Degrazia, G.A., Tirabassi, T., Rizza, U., 2002. Evaluation of a new eddy diffusivity parameterization from turbulent Eulerian spectra in different stability conditions. Atmos. Environ. 36, 67–76. Moreira, D.M., Moret, M., 2018. A new direction in the atmospheric pollutant dispersion inside of the Planetary Boundary Layer. J. Appl. Meteorol. Climatol. 57 (1), 185–192. Moreira, D.M., Santos, C.A.G., 2019. New approach to handle gas-particle transformation in air pollution modelling using fractional derivatives. Atmos. Pollut. Res. 10, 1577–1587. Moreira, D.M., Vilhena, M.T., 2009. Air Pollution and Turbulence: Modeling and Applications. CRC Press, Boca Raton, Florida, pp. 354. Moreira, D.M., Vilhena, M.T., Buske, D., Tirabassi, T., 2009. The state-of-art of the GILTT method to simulate pollutant dispersion in the atmosphere. Atmos. Res. 92, 1–17. Moreira, D.M., Tirabassi, T., Carvalho, J.C., 2005c. Plume dispersion simulation in low wind conditions in stable and convective boundary layers. Atmos. Environ. 39 (20), 3643–3650. Moreira, D.M., Tirabassi, T., Vilhena, M.T., Goulart, A.G., 2010. A multi-layer model for pollutant dispersion with dry deposition to the ground. Atmos. Environ. 44 (15), 1859–1865. Moreira, D.M., Moraes, A.C., Goulart, A.G., Albuquerque, T.T., 2014. A contribution to solve the atmospheric diffusion equation with eddy diffusivity depending on source distance. Atmos. Environ. 83, 254–259. Moreira, D.M., Vilhena, M.T., Carvalho, J.C., Degrazia, G.A., 2005a. Analytical solution of the advection-diffusion equation with nonlocal closure of the turbulent diffusion. Environ. Model. Softw 20 (10), 1347–1351. Moreira, D.M., Vilhena, M.T., Tirabassi, T., Buske, D., Cotta, R.M., 2005b. Near source atmospheric pollutant dispersion using the new GILTT method. Atmos. Environ. 39,

10