In this study, an analytical solution for the steady-state fractional advection–diffusion equation was obtained to simulate the atmospheric dispersion of pollutants in a vertically inhomogeneous planetary boundary layer. The authors propose a method that uses the modified generalized integral Laplace transform technique to solve the transformed problem with a fractional derivative, resulting in a more general solution. The model results were compared with the fractional Gaussian model and demonstrate that, when considering an experimental dataset under moderately unstable conditions, fractional-derivative models perform better than traditional integer-order models.
There currently is strong interest in analytical solutions of differential equations, in which the fractional calculation, especially fractional differential equations, has been used to describe many natural phenomena (Zaslavsky 1994; Meerschaert and Tadjeran 2004; Gorenflo and Mainardi 2009; Schumer et al. 2009). This calculation has only recently obtained practical applications, however. In this way, the fractional calculation became a very useful tool to study anomalous diffusion and other transportation processes, initially with more attention to the hydrological and porous environment (Benson et al. 2013; Deseri and Zingales 2015). Debnath (2003) contains a good description of the recent applications of fractional calculus to science and engineering. Furthermore, we cite the books of Miller and Ross (1993) and Podlubny (1998).
The transport of a pollutant under the combined effects of advection and diffusion is described by the diffusion–advection equation, with which many physical processes can be modeled. Atmospheric pollution caused by natural or anthropogenic effects is an example that has been systematically modeled by traditional differential equations (integer order) (Lin and Hildemann 1997; Moreira et al. 2005a,b, 2009, 2014; Sharan and Modani 2006; Essa et al. 2007; Tirabassi et al. 2008; Perez-Guerrero et al. 2012; Pimentel et al. 2014). There is, however, a gap in the analytical solutions of differential equations of fractional order in the literature regarding atmospheric pollutant dispersion, where the use of fractional calculus in the modeling of turbulent diffusion is justified by the nondifferentiable behavior in the problem and by the presence of anomalous diffusion. In this way, the literature using fractional calculus to model atmospheric diffusion problems is practically nil. Furthermore, the literature presents a trend in the use of numerical solutions of fractional differential equations in several problems, among them the Navier–Stokes equation (Birajdar 2014; Guo et al. 2015). Analytical solutions of fractional differential equations applied to practical problems still require more attention, however. Although this work proposes the extension of a well-known method (Buske 2004; Moreira et al. 2015) with the inclusion of fractional derivatives in the advection–diffusion equation, the main idea is to build a road to an analytical solution generalized to the fractional advection–diffusion equation in atmospheric problems.
In this context, Goulart et al. (2017) proposed an alternative in the understanding of the atmospheric dispersion process by modifying the mathematical structure of the advection–diffusion equation, introducing fractional operators in the equation that governs the distribution of contaminants in the atmosphere. This modification more realistically represents the spatial evolution of the concentration of atmospheric pollutants dispersed in a turbulent flow, resulting in an equation with a fractional spatial derivative in the advection term, given by
where is the crosswind-integrated concentration, C represents the Caputo derivative, α represents the order of the fractional spatial operator, u is the longitudinal mean speed, and KZ is the vertical eddy diffusivity. Equation (1) is the general steady-state fractional advection–diffusion equation displaying the anomalous power-law mean-squared displacement (Goulart et al. 2017).
At this point it is important to mention that the phenomenon of turbulent diffusion is one of the challenges of modern science, for which the mathematical model represented by traditional integer-order differential equations does not adequately describe the problem of turbulent diffusion. Thus, it is expected that the usual advection–diffusion equation does not represent the anomalous dispersion of pollutants so that the usual procedure in the literature to modify the Eulerian models is to assume that the physical structure of the turbulent flow and velocity field are described by eddy diffusivities and average velocities, which are considered to be functions of spatial coordinates. That is, in the traditional models of pollutant dispersion, the parameters u and KZ are usually considered to be functions of x and z variables so as to handle the anomalous diffusion introduced by the turbulence. In particular, the memory effect, important in near regions from high and low continuous sources, has been introduced in many works through the vertical eddy diffusivity (Degrazia et al. 2001; Moreira et al. 2014). Unlike the usual case, here the memory effect has been introduced through the fractional derivative. In this work, because of the modification of the mathematical structure of the traditional advection–diffusion equation (Goulart et al. 2017), we propose an advance in the understanding of the dispersion process of atmospheric pollutants with the obtaining of an analytical solution coming from the generalization of the method of the generalized integral Laplace transform technique (GILTT), which will now be represented with the solution of the fractional advection–diffusion equation in atmospheric problems. Despite the apparent “small modification” in the equation with the introduction of the fractional derivative represented by the parameter α, this term represents an important step in the modeling of the anomalous dispersion process of atmospheric pollutants. Following the objective of the work, there is currently a need to improve or adapt the methods and techniques of existing analytical solutions from the classical case to the fractional one. At first glance, this process seems simple and straightforward, but this process is complicated and requires attention because the fractional calculation needs additional conditions to be defined correctly.
Therefore, the novelty of the study presented here is to solve Eq. (1) analytically, with u and KZ depending on the vertical z variable, using the GILTT (Wortmann et al. 2005; Moreira et al. 2005a, 2009), in which the transformed problem arising from this method contains the fractional derivative, known in this study as the method of the α-GILTT. This procedure is used for the first time in the field of atmospheric pollutant dispersion. In addition, different from the solution obtained in the work of Goulart et al. (2017), this method is general. In this context, it is emphasized that, for more complex situations of dispersion of pollutants in the atmosphere, we must resort to numerical models, but the analytical solutions are always necessary for a quick analysis; besides easily handling the primary variables involved in the process of modeling, they serve as an excellent tool for the evaluation of the results from the most sophisticated numerical models. Furthermore, we can say that, with this procedure, the analytical solution of the fractional advection–diffusion equation obtained here represents a powerful tool to calculate the concentration field of pollutants, enabling a better understanding of the dispersion process with less computational effort because of the analytical characteristics of the model (Arya 1999; Moreira and Vilhena 2009).
2. Fractional-approach solution
In this work, we consider that u and Kz depend on variable z; therefore, we may write Eq. (1) as follows:
Thus, to solve the problem represented in Eq. (2) with the GILTT, the equation can be rewritten as follows:
with the following boundary conditions:
where h is the height of the planetary boundary layer (PBL), z0 is the roughness length, Q is the emission rate, Hs is the source height, and δ() is the Dirac delta function.
Following the method of the GILTT, we can now construct the problem associated with the Sturm–Liouville theory:
and this problem has the traditional solution given by
where the eigenvalues λi and the eigenfunctions Ψi(z) satisfy the orthonormality condition:
with Nm given by
Therefore, assuming the above considerations, we can write the inverse GILTT formula with the following solution:
Thus, by inserting this equation into Eq. (3), we obtain
It should be remembered that we have adopted a prime to indicate the first derivative and a double prime for the second derivative and that α is the order of a fractional derivative. Following the GILTT formalism, we multiply Eq. (8) by and integrate over z from 0 to h:
for j = 0, 1, 2, …, N. For convenience, it is possible rewrite Eq. (9) in matrix form, resulting in
where Y(x) is a vector, the matrix is defined as = −1, and α is an arbitrary real number. At this point, it is important to note that the novelty of this work appears here: the transformed problem with the fractional derivative.
where and are constant matrices of dimension N × N and Y0 is a known vector with N components.
Thus, to solve the problem represented by Eq. (10), following the GILTT formalism, we apply the Laplace transform technique on variable x. We recall that the Laplace transform (indicated by curly braces) of a fractional derivative of α order is given by Caputo’s formula, as follows:
where s is the transformed variable. Therefore, in our case, we have the following:
Hence, Eq. (10) can be written as
where results from the application of the Laplace transform, and the matrix can be written as
where is the identity matrix. Thus, Eq. (15) becomes the following:
Applying the inversion of Laplace transform on Eq. (18), we have
where L−1 is the inverse Laplace transform. Next, using the standard theory for Laplace transform, we obtain
where di are the eigenvalues of matrix , and Eα is the Mittag-Leffler function given by
where Γ() is the gamma function. The Mittag-Leffler function is intrinsic to the solution of equations with fractional derivatives. Last, substituting Eq. (20) into Eq. (19), we obtain the following solution to the problem given by Eq. (19):
We can readily observe that the solution to the transformed problem given by Eq. (22) is more generic, in the sense that it returns the traditional GILTT solution when assuming α = 1. Recall that the elements of matrices and are respectively given by
resulting in the following source condition:
For the vertical eddy diffusivity, the formulation of Troen and Mahrt, described in the work of Pleim and Chang (1992), is used throughout the whole PBL and is defined as follows:
where h is the height of the PBL, w* is the convective velocity, and k is the von Kármán constant and is approximately 0.4. The convective velocity from the experimental dataset that is used in this work is provided from the expression w* = u*[−h/(kL)]1/3 (Degrazia et al. 2001), where L is the Monin–Obukhov length.
The model calculates the mean wind velocity using the following similarity equations (Seinfeld and Pandis 1998):
where z0 is the roughness length of the terrain, zb = min(|L|, 0.1h), and Ψm is the stability function expressed by Paulson (1970).
3. Numerical results
To verify the model’s physical consistency, this paper takes as its basis the traditional experimental data from Copenhagen that are described in the work of Gryning and Lyck (2002). This experiment is considered to represent moderate convection (Pasquill and Smith 1983).
We compare the results obtained by the fractional GILTT model (here called α-GILTT) and the fractional Gaussian model (GM) from the work of Goulart et al. (2017) (here called α-GM). In the work of Goulart et al. (2017), a diffusion coefficient that is dependent on the longitudinal distance was used, which was taken as an average value in this direction. Since the Gaussian model is for constant coefficients, we adopted the diffusion coefficient represented in this work by Eq. (28) as the mean value in the vertical direction (height of the PBL).
The micrometeorological data used in the initial simulations of this study are those from experiment 9 at Copenhagen, with u* = 0.75 m s−1 (friction velocity), w* = 1.9 m s−1 (convective velocity), h = 2090 m (height of the PBL), and L = −289 m (Monin–Obukov length). These micrometeorological parameters can be obtained from the work of Degrazia et al. (2001).
Figure 1 shows the graph of the concentration integrated laterally, at ground level, as a function of the distance from the source using data from experiment 9 at Copenhagen. We observe that the concentration generated by the α-GILTT model with α = 0.95 provides the best results with the parameterization used. The concentration for the distance of 2100 m was best represented using α = 0.9. Parameter α does not change the value of the peak concentration, which is one of the most important aspects in the context of air pollution, but it does modify its position. In comparing this case with the α-GM, we observe that, for the range of tracer collection points in this experiment (2–6 km), the results are similar but the results from the fractional GILTT model are still better. Indeed, the α-GM shows a much higher peak than the fractional GILTT at short distances, where the concentration’s maximum position practically does not change.
In its sequence, Fig. 2 shows the vertical concentration profile for different distances from the source: 0.25, 1, 10, and 20 km. This qualitative analysis is used to detect possible differences between the models. We readily observe the influence of the α term on both models in Fig. 2. For the distance closest to the source (0.25 km), a concentration peak is observed in the region of the source’s height that is stronger in the α-GILTT model and intensifies when parameter α decreases. The concentration at ground level generated by the α-GM is higher (approximately 4 times as great with α = 1) for this distance. As the distance from the source increases, we observe a reduction of the peak and of the concentration’s numerical value with a tendency for vertical homogenization. This behavior is most clearly observed for higher distances (10 and 20 km). At the distance of 20 km, complete homogenization of the concentration in the vertical direction is observed, both for α = 1 and for α = 0.95. For α = 0.85, complete homogenization in the vertical direction has not been observed to date. At all simulated distances that we observe, there is a clear distinction between the results of the α-GM and the results of the α-GILTT model, which considers continuous profiles of the wind velocity and diffusion coefficient parameters. Numerical convergence of the proposed solution for the concentration with increasing numbers of eigenvalues is reached considering the traditional integer atmospheric diffusion equation and the fractional diffusion equation for different α parameters. From a practical point of view, we used N = 100.
For a more general analysis comparing the models, we use statistical indices that are well known in the literature (Hanna 1989). Table 1 shows the results of the statistical indicators (see the table caption for their definitions) obtained from the α-GM and the α-GILTT model, where different values for parameter α have been considered.
Table 1 shows that the results obtained with the inclusion of the parameterization dependent on the variable z [u(z) and K(z)] in the proposed dispersion model (α-GILTT), with α = 0.96 (in boldface font), showed a strong correlation with the experimental data because the results included a correlation coefficient of 0.9, Fat2 of 100%, and the lowest NMSE (0.05). The α-GM proposed by Goulart et al. (2017), who considered constant values for u and K, also showed satisfactory results, with the best indices for α = 0.92 (in boldface font). In general, the statistical results are comparable for almost all simulations with α > 0.90 because they intersect at the confidence limit (not shown), except when α = 0.85, for both models. It is noted that for the Copenhagen experiment, which has experimental data for concentration and distances between 2 and 6 km, near the ground (z ≈ 0) the best results are similar, with a slight advantage for the α-GILTT solution. In the figures, we observe that the concentration is strongly dependent on the vertical component; in addition, the models show very different results for distances of less than 1.5 km and greater than 6 km, where in the α-GILTT model the concentration decreases more rapidly with increasing distance from the source.
In this work, we proposed and investigated the sensitivity of an analytical solution of the two-dimensional steady-state fractional advection–diffusion equation for estimating the concentrations of pollutants in the PBL. The model allows the use of realistic wind profiles and vertical eddy diffusivity in any atmospheric stability condition, although this work reports on simulations in moderately unstable conditions.
The simulations and comparisons performed with a Gaussian dispersion model suggest better statistical results for the α-GILTT model when using parameter α with this work’s proposed parameterization. The performance of the Gaussian model does present comparable results with α = 0.92 for the Copenhagen experimental data that are only for ground-level concentrations. This statement is demonstrated by the statistical indicators. It is important to mention that the simulations are very different for greater heights, however.
The obtained solution is analytical, in the sense that no approximations have been made in its derivation, apart from the fact that the transformed problem is solved by the analytical inversion of the Laplace transform, resulting in a more general solution with the introduction of the Mittag-Leffler function, which is intrinsically associated with fractional equations. The solution of the fractional problem gives the traditional solution (integer order) when α = 1. We hope that with the formulation reported we have paved the way to update the GILTT formulation to simulate pollutant dispersion considering anomalous dispersion, which is an important advance in the atmospheric sciences. We shall focus our future attention in this direction.
The authors thank the National Council for Scientific and Technological Development (CNPq) and the Bahia State Research Support Foundation (FAPESB).