Generalized Skewed Model for Spatial-Fractional Advective–Dispersive Phenomena

Abstract

The conventional mathematical model expressed by the advection–dispersion equation has been widely used to describe contaminant transport in porous media. However, studies have shown that it fails to simulate early arrival of contaminant, long tailing breakthrough curves and presents a physical scale-dependency of the dispersion coefficient. Recently, advances in fractional calculus allowed the introduction of fractional order derivatives to model several engineering and physical phenomena, including the anomalous dispersion of solute particles. This approach gives birth to the fractional advection–dispersion equation. This work presents new solutions to the fractional transport equation that satisfies the initial condition of constant solute injection in a semi-infinite medium. The new solution is derived based on a similarity approach. Moreover, laboratory column tests were performed in a Brazilian lateritic soil to validate the new solution with experimental data and compare its accuracy with the conventional model and other fractional solutions. The new solution outperforms the existing ones and reveals an interesting fractal-like scaling rule for the diffusivity coefficients.

1. Introduction

Understanding how solute transport occurs inside a porous medium is one of the most important problems to geo-environmental engineers. Natural porous media exhibit heterogeneity at all scales, mainly due to the geologic formation processes to which they have been subjected. These processes involve fractures, faults and cracks having various structures, discontinuities, slopes, propagations, shapes and geometries, arrangements and sizes [1,2].

The widely known mathematical model used to describe contaminant transport phenomenon on porous media is expressed by the advection–dispersion equation (ADE):

$$\frac{\partial C}{\partial t}=-v\frac{\partial C}{\partial x}+D\frac{{\partial }^{2}C}{\partial x^2}$$

(1)

where C is the solute concentration $\left[M{L}^{-3}\right]$, v the flow velocity $\left[L{T}^{-1}\right]$ and D the dispersion coefficient $\left[L^2{T}^{-1}\right]$, t is time $\left[T\right]$ and x stands for the single spatial variable of the problem $\left[L\right]$.

Equation (1) has been extensively studied in literature and solved for several relevant initial and boundary conditions to simulate different scenarios [3]. For example, Ogata and Banks [4] describe the modelling of a continuous source of contamination injected at the top of a soil layer by solving the ADE for a semi-infinite medium instead of a $x=-\infty$ to $x=\infty$ boundary condition.

Even though the ADE can be used to model the solute transport process in homogeneous porous media, it fails to capture some important features present in heterogeneous porous media [5]. Studies indicate that such inadequacy comes from the fact that the ADE is based on assuming the applicability of Fick’s first law, while the solute transport process in the heterogeneous porous media follows a non-Fickian or anomalous mechanism [5]. In general, Fick’s first law indicates how the diffusive flux can be related to the gradient of the concentration. More specifically, it indicates that the flux goes from high concentration regions to low concentration regions, with a magnitude that is proportional to the concentration gradient (partial derivative with respect to space). Put in a simple manner, this law states that a solute will move from a region of high concentration to a region of low concentration across a concentration gradient [6].

On the other hand, non-Fickian transport gives place to some interesting behaviours, e.g., early arrival of contaminant, long tailing breakthrough curves (BTCs), and scale-independency of the dispersion coefficient are among the observed non-Fickian or anomalous characteristics of real-world contaminant transport [7,8].

More recently, advances in fractional calculus theory allowed the introduction of fractional derivatives in a series of differential equations used to simulate various phenomena, including contaminant transport in porous media. Different approaches to a fractional transport equation are available in the literature, such as the inclusion of a fractional definition in the time derivative [9], in the time and space derivatives [10] and in both the advection and dispersion space derivatives [11].

By considering the inclusion of fractional derivatives only with respect to space, Benson et al. [12] proposed the following fractional advection–dispersion equation (FADE):

$$\frac{\partial C}{\partial t}=-v\frac{\partial C}{\partial x}+\frac{1}{2}(1+\beta )\mathcal{D}\frac{{\partial }^{\alpha }C}{\partial x^{\alpha }}+\frac{1}{2}(1-\beta )\mathcal{D}\frac{{\partial }^{\alpha }C}{\partial {(-x)}^{\alpha }}$$

(2)

where $\alpha$ is the derivation index, $\beta$ the asymmetry parameter, ${\partial }^{\alpha }/\partial x^{\alpha }$ is the forward Riemann–Liouville fractional derivative [13,14], ${\partial }^{\alpha }/\partial {(-x)}^{\alpha }$ is the backward Riemann–Liouville fractional derivative [15] and $\mathcal{D}$ the fractional dispersion coefficient $\left[L^{\alpha }{T}^{-1}\right]$.

The FADE is a spatially non-local model which is capable of describing movements of solute particles that experience very large transitions [16]. These large transitions may be due to high heterogeneity, very long spatial autocorrelation, and high-velocity preferential flow paths, resulting in a long leading tail and, consequently, a super-diffusive transport process [12,16].

Using space-fractional approaches gives rise to stable distributions for concentration fields, which jeopardizes an interpretation in terms of statistical moments. In this particular case, the mean squared displacement diverges, which is a nonphysical property. This issue may be remedied by working with applying cutoffs to the integrals involved and by considering only a finite number of sampled tracer particles in a finite range, during a finite time window. This leads to a truncated Lévy distribution with finite moments, which allows this type of model to be validated in physical terms as well [17].

The literature reveals the proposition of other generalizations of the transport models [18]. For example, generalizations with long-tailed memory due to broad distributions in trapping times have been discussed [17,19]. Indeed, as can be seen in [20], there exists direct, constructive evidence for long-tailed transition time distributions. Further studies in that direction can also be seen in [21,22]. Anyway, we focus our studies on Equation (2) and its solutions as it has been extensively used in soil sciences.

This way, Equation (2) has been solved by Benson et al. [12] for the step function initial condition as well as for a constant source condition. Besides, it has been applied to model real experimental data [7,23,24,25].

Fractional derivatives are non-local, meaning the fractional derivative of a function at a given point depends on the characteristics of the function across the entire spatial domain [26]. Therefore, as mentioned above, the FADE is spatially non-local, signifying that the concentration change at some location in time depends not only on the current position but also on the properties of a wide range of upstream and downstream locations [16]. Theoretically, the diffusivity coefficient of the FADE is scale-independent and the derivative index $\alpha$ reflects the scale effects [27]. Besides, Mehdinejadiani et al. [28] indicated that there exists a relation between the fractional differentiation order $\alpha$ of a fractional model and the heterogeneity degree of a given porous medium.

Even though in theory the diffusivity coefficient of the FADE is scale-independent, the literature reveals that experiments indicate that such parameter varies with scale. Lu et al. [29] considered both the FADE and the ADE to model solute transport in a highly heterogeneous aquifer. Their results indicate that despite the FADE being superior to the ADE (in terms of fitting the data), both the dispersion coefficients of the FADE and of the ADE presented scale dependency.

The same tendency was verified by Huang et al. [30]. These authors showed that the application of the FADE to model solute transport in saturated porous media did not eliminate the scale effect on the diffusivity coefficient, but significantly reduced it. More recently, Moradi and Mehdinejadiani [16] performed a massive experimental campaign to study the scale dependency of the diffusivity coefficient, which verified the same characteristics.

All the works previously discussed considered the symmetrical version of Equation (2), i.e., when $\beta =0$. This comes from the fact that the analytical solutions presented in the literature consider such simplifying assumptions [7,12].

In the present paper, a novel skewed (asymmetrical) solution is proposed for Equation (2) with $-1\le \beta \le 1$. The new solution is presented to improve the modelling of a continuous source of contamination injected at the top of a soil layer. Therefore, a semi-infinite boundary condition is considered. Additionally, laboratory column tests were conducted to verify the behaviour and precision of the solution hereby presented. Finally, the scale-dependency of the diffusivity coefficient obtained for the tests performed is discussed.

2. Known Solutions to the ADE and to the FADE

At first, one shall consider the solution to the 1D ADE, expressed in Equation (1), as given by Ogata and Banks [4] in its full form:

$$\frac{C}{C_0}=\frac{1}{2}\left[\mathrm{erfc}\left(\frac{x-vt}{2\sqrt{Dt}}\right)+e^{\frac{vx}{D}}\mathrm{erfc}\left(\frac{x+vt}{2\sqrt{Dt}}\right)\right]$$

(3)

which satisfy the semi-infinite medium boundary conditions:

$$\frac{C}{C_0}=\left\{\begin{array}{c}1,x=0;t>0\hfill \\ 0,\forall x;t=0\hfill \\ 0,x=\infty ;t\ge 0\hfill \end{array}\right.$$

(4)

In general, solutions to common solute transport boundary value problems (BVP) are obtained through Laplace or Fourier transforms. This was the approach considered by Ogata and Banks [4] and Benson et al. [12].

In fact, the continuous source BVP using the classical ADE is reasonably approximated by Ogata and Banks [4]:

$$\frac{C}{C_0}=\frac{1}{2}\left[1-\mathrm{erf}\left(\frac{x-vt}{2\sqrt{Dt}}\right)\right]$$

(5)

Ogata and Banks [4] discuss that the second term in (3) may be disconsidered for small values of D or big values of distance x, as solution (3) becomes equal to (5) and the boundary condition equivalent to an infinite medium $C(-\infty ,t)=C_0$. They indicate that for $x>500D/v$, the second term can be removed with a maximum error of about 3%. Figure 1 shows the effect of considering an infinite medium boundary condition, which creates an always symmetrical solution (around $C/C_0=0.5$), and a semi-infinite medium boundary condition that becomes asymmetrical as $x\to 0$ to correctly simulate the initial concentration $C_0$ at the top of a soil layer.

By considering a similarity approach, Benson et al. [12] analysed (5) and proposed the following formula for the solution of the corresponding FADE:

$$\frac{C}{C_0}=\frac{1}{2}\left[1-{\mathrm{serf}}_{\alpha }\left(\frac{x-vt}{{\left(\mathcal{B}t\right)}^{1/\alpha }}\right)\right]$$

(6)

where $\mathcal{B}=|cos(\pi \alpha /2)|\mathcal{D}$ and ${\mathrm{serf}}_{\alpha }$ is the $\alpha$-stable error function, mathematically defined as twice the integral of a symmetric $\alpha$-stable density,$f\left(x\right|\alpha ,0,0,1)$, from 0 to the argument z, i.e.,:

$${\mathrm{serf}}_{\alpha }\left(z\right)=2{\int }_0^{z}f\left(x\right|\alpha ,0,0,1)dx$$

(7)

Schneider [31] gave a representation of the probability density function (PDF), $f\left(x\right|\alpha ,\beta ,\delta ,\gamma )$, of a $\alpha$-stable Levy random variable in terms of the H-function, which is a generalized hypergeometric special function defined in terms of a contour integral of products and quotients of Gamma functions [32]. For simplicity, if a random variable X follows an $\alpha$-stable distribution with parameters $\alpha$, $\beta$, $\delta$ and $\gamma$, we say $X\sim S(\alpha ,\beta ,\delta ,\gamma )$.

Expressions equivalent to the ones obtained by Schneider [31] for the PDF of $X\sim S(\alpha ,\beta ,\delta ,\gamma )$ with $0<\alpha <2$, $\alpha \ne 1$, were presented by Rathie et al. [33]. Additionally, Rathie and Ozelim [34] provided expressions for the cumulative distribution function (CDF), $F\left(x\right|\alpha ,\beta ,\delta ,\gamma )$, of an $\alpha$-stable random variable in terms of the H-function. Since the numerical evaluation of the H-function became available only recently, computational software such as Mathematica [35] have other built-in implementations for both the PDFs and CDFs of $\alpha$-stable random variable. In fact, in the present paper, Mathematica shall be used to evaluate the equations hereby presented.

Equation (6) may be rearranged in terms of the CDF of an $\alpha$-stable standard random variable by noticing that, when $\beta =0$, $\delta =0$ and $\gamma =1$:

$$\begin{array}{ccc}\hfill F\left(x\right|\alpha ,0,0,1)& =& {\int }_{-\infty }^{x}f\left(x\right|\alpha ,0,0,1)dx\hfill \\ & =& {\int }_{-\infty }^{0}f\left(x\right|\alpha ,0,0,1)dx+{\int }_0^{x}f\left(x\right|\alpha ,0,0,1)dx\hfill \\ & =& \frac{1}{2}+\frac{{\mathrm{serf}}_{\alpha }\left(x\right)}{2}\hfill \end{array}$$

(8)

Thus, from Equation (6) and (8):

$$\frac{C}{C_0}=1-F\left(\left(\frac{x-vt}{{\left(\mathcal{B}t\right)}^{1/\alpha }}\right)|\alpha ,0,0,1\right)$$

(9)

The representation in (9) was also pointed out in Benson et al. [12]. Due to the scaling and shifting properties of stable random variables, Equation (9) may be rewritten as:

$$\frac{C}{C_0}=1-F\left({x|\alpha ,0,vt,\left(\mathcal{B}t\right)}^{1/\alpha }\right)$$

(10)

It is clear from (9) that skewness is not present in the above model. By further analysing Equation (9) a new solution to the FADE is proposed. It is also important to consider other solutions presented in the literature. Therefore,

Table 1. Previous models reported in the literature.

References	Type of Solution
Benson et al. [36]	Analytical, with data application/comparison, $\beta =0$
Meerschaert et al. [37]	Analytical, continuous time random walk (CTRW) approach, $\beta =0$
Lu et al. [29]	Analytical, 3D space fractional derivatives, $\beta =1$
Cortis et al. [38]	Numerical, CTRW approach, $\beta$ not assignable directly
Deng et al. [39]	Numerical, space fractional derivative, $\beta =0$
Tadjeran et al. [40]	Numerical, with space fractional derivative, for $v=0$ and $\beta =0$
Burrage [41]	Numerical, with fractional derivatives in space and time, for $v=0$ and $\beta$ not assignable directly
Fomin et al. [42]	Numerical, with data application/comparison, $\beta =0$
Ouloin et al. [43]	Numerical, space fractional derivative, $\beta =0$
Saffarian and Mohebbi [44]	Numerical, 2D space fractional derivatives, $\beta$ not assignable directly

has been created to indicate that none of the previous models reported in the literature provide a simple and easy-to-use expression when asymmetry is present.

3. New Solution to the ADE and to the FADE

When deriving (6), Benson et al. [12] used a similarity approach. In that case, they introduced a $\alpha$-stable error function which would reduce to the Gaussian case when $\alpha =2$. In the present paper, a similar approach shall be used to obtain a new solution to the ADE and the FADE.

3.1. ADE Approximate Similarity Solution

Equation (3) fully satisfies the boundary conditions imposed to the PDF. In special, it is clear from (3) that for $x=0$, $C=C_0$. This comes from the fact that the value of the functions inside the square brackets in the right-hand side of (3) is 2. This value is then normalized by the $1/2$ factor outside the brackets, resulting in $C/C_0=1$ when $x=0$. By noticing that, (3) may be rewritten as:

$$\frac{C}{C_0}=\frac{\left[\mathrm{erfc}\left(\frac{x-vt}{2\sqrt{Dt}}\right)+e^{\frac{vx}{\mathcal{D}}}\mathrm{erfc}\left(\frac{x+vt}{2\sqrt{Dt}}\right)\right]}{{\left[\mathrm{erfc}\left(\frac{x-vt}{2\sqrt{Dt}}\right)+e^{\frac{vx}{D}}\mathrm{erfc}\left(\frac{x+vt}{2\sqrt{Dt}}\right)\right]}_{x=0}}$$

(11)

Equations of the type of (11) will always satisfy the boundary conditions in (4) when the function in the numerator is monotonically decreasing with respect to x.

One important consequence of disconsidering the second term in (3) is that this approximate solution (5) no longer satisfies the boundary conditions in (4).

In fact, at $x=0$, (5) provides:

$$\frac{C}{C_0}=\frac{1}{2}\left[1-\mathrm{erf}\left(\frac{-vt}{2\sqrt{Dt}}\right)\right]$$

(12)

This way, by considering the procedure given in (11), a new approximate solution to the ADE would be of the form:

$$\frac{C}{C_0}=\frac{\mathrm{erfc}\left(\frac{x-vt}{2\sqrt{Dt}}\right)}{\mathrm{erfc}\left(\frac{-vt}{2\sqrt{Dt}}\right)}$$

(13)

The solution in (13) satisfies the boundary and initial conditions of the original solution in (4).

3.2. FADE Approximate Similarity Solution

On the other hand, in order to build a new approximate solution to the FADE, one shall consider (10). By applying the same normalizing concept as applied in (11), the new approximate solution to the FADE may be written as:

$$\frac{C}{C_0}=\frac{1-F\left({x|\alpha ,\beta ,vt,\left(\mathcal{B}t\right)}^{1/\alpha }\right)}{1-F\left({0|\alpha ,\beta ,vt,\left(\mathcal{B}t\right)}^{1/\alpha }\right)}$$

(14)

Figure 2 shows that by considering the normalization factor different than 2, the curve is pulled up to meet the boundary conditions at $x=0$. The differences between both solutions are, as expected, more relevant when x tends to 0. Additionally, the possibility to consider asymmetric diffusion by using $\beta \ne 0$ provides a new perspective on the modelling of contaminant transport in porous media.

It is worth noting that (14) is an approximate solution only for the advection–dispersion phenomena. For pure diffusion, i.e., when $v=0$, solution (14) is exact and equivalent to (3) for $\alpha =2$.

Furthermore, it is possible to remove the factor $|cos(\pi \alpha /2)|$ from the fractional dispersion coefficient $\mathcal{B}$ by using the fractional derivative definition based on the Riesz potential instead of the Riemann–Liouville case. The Riesz potential derivative can be given as [45,46]:

$$\frac{{\partial }^{\alpha }g\left(x\right)}{\partial x^{\alpha }}=-\frac{\frac{{\partial }^{\alpha }g\left(x\right)}{\partial x^{\alpha }}+\frac{{\partial }^{\alpha }g\left(x\right)}{\partial {(-x)}^{\alpha }}}{2cos(\alpha \pi /2)}$$

(15)

Thus, the resulting Riesz-fractional advection-dispersion equation (R-FADE) is given as:

$$\frac{\partial C}{\partial t}=-v\frac{\partial C}{\partial x}-\frac{(1+\beta )}{2cos\left(\pi \alpha /2\right)}\mathcal{D}\frac{{\partial }^{\alpha }C}{\partial x^{\alpha }}-\frac{(1-\beta )}{2cos\left(\pi \alpha /2\right)}\mathcal{D}\frac{{\partial }^{\alpha }C}{\partial {(-x)}^{\alpha }}$$

(16)

For the Dirac delta BVP, Equation (16) is solved using Fourier transforms (Appendix A) similarly to the FADE solution in Benson et al. [12]. In the same manner as (14), the approximate solution of the R-FADE that satisfies the initial and boundary conditions in (4) is given by:

$$\frac{C}{C_0}=\frac{1-F\left({x|\alpha ,\beta ,vt,\left(\mathcal{D}t\right)}^{1/\alpha }\right)}{1-F\left({0|\alpha ,\beta ,vt,\left(\mathcal{D}t\right)}^{1/\alpha }\right)}$$

(17)

Figure 3 shows the effect of factor $|cos(\pi \alpha /2)|$ on solutions (14) and (17), specifically the behaviour of the scale parameter γ. Considering the behaviour of a fractional order root, when its argument is greater than 1 (Figure 3a,b), a decrease in fractional index $\alpha$ causes an increase in scale parameter $\gamma$ and thus represents faster particle dispersion. As $\alpha$ tends to 1, the factor $|cos(\pi \alpha /2)|$ in (14) zeroes the scale parameter, only presenting the inversely proportional relation to $\alpha$ for bigger values of $\mathcal{D}t$ or $\alpha$ values closer to 2. On the other hand, when the fractional root argument $\mathcal{D}t$ is smaller than 1 (Figure 3c,d), $\gamma$ becomes proportional to $\alpha$ and solution (17) describes a gradual decreasing function as $\alpha \to 1$ while (14) quickly goes to zero.

Since (17) is a semi-analytical solution, i.e., it needs numerical values to calculate the CDF of the $\alpha$-stable distribution, a Python 3 [47] script is provided in Appendix B to facilitate the understanding and application of the new solution.

Now that the new solutions have been obtained, one shall discuss the advantages of including the skewness parameter into the solution.

3.3. Skewness Parameter

The addition of the asymmetry parameter $\beta \ne 0$ in solutions (14) and (17) allows one to simulate preferential anomalous dispersion in one direction or another, resulting in asymmetric breakthrough curves (BTCs) that may better represent experimental data.

For $\beta =1$, i.e., maximum positive skewness, the backward fractional derivative in (2) and (16) is multiplied by zero and only the forward fractional derivative is considered. On the other hand, for $\beta =-1$, i.e., maximum negative skewness, the forward fractional derivative is multiplied by zero, resulting in only the backward fractional derivative. Figure 4 shows the new solution for a range of $\beta$ values.

To better understand the effect of positive and negative skewness in the fractional dispersion phenomenon, Schumer et al. [48] point out the convergence of a random walk with $\alpha$-stable random variable $X_{\alpha }$ to the chaotic process of Lévy flight and its equivalence to the fractional advection–dispersion solution. As such, Figure 5 and Figure 6 simulate the walk of a random particle according to an $\alpha$-stable distribution with fractional index $\alpha =1.2$ and asymmetry parameter $\beta =1$ (Figure 5) and $\beta =-1$ (Figure 6). The fractional index $\alpha$ indicates the probability of larger jumps at each step taken by the particle, where values closer to 1 indicate a higher frequency of anomalous movement, resulting in a faster spread of particles (super-diffusion). The asymmetry parameter $\beta$ indicates the preferential direction of such movement, where a positively skewed distribution favours anomalous steps in the positive direction, and a negatively skewed distribution favours anomalous steps in the negative direction, resulting in an asymmetrical spread of particles.

This way, asymmetry may enhance the modelling capability of the equations proposed while compared to standard symmetrical models. In the next section, experimental data are fitted to this new solution to show its applicability.

4. Application of the New Solution to Experimental Data

To compare the efficiency and precision of the new solution with respect to the ones presented by Benson et al. [12] and Ogata and Banks [4], laboratory column tests were performed and their respective breakthrough curves fitted with each model.

The soil used in the tests is a lateritic inorganic sandy clay with $n=56\%$ porosity, collected in Minas Gerais, Brazil. The soil was preserved in its natural conditions, in order to preserve the original void disposition, and moulded into cylinders with heights of 6 cm and 12 cm and diameters of 5 cm. These cylinders (columns) were then percolated with a multi-species solution containing copper (${\mathrm{Cu}}^{2+}$) and lead (${\mathrm{Pb}}^{2+}$) chloride [49]. A total of four 6 cm samples were tested with a ${\mathrm{Cu}}^{2+}$ contaminant solution, while two samples of the same height were subjected to a ${\mathrm{Pb}}^{2+}$ contaminant solution. To assess the scale effects (how the new diffusion equation can account for flow paths with different lengths and flow networks), four samples with heights of 12 cm were tested considering a ${\mathrm{Cu}}^{2+}$ contaminant solution. Figure 7 presents a detail of the soil cylinders (columns), while the general experimental setup can be seen on Figure 8.

The resulting breakthrough curves were fitted to each solution using a non-linear least-squares regression method (LSM):

$$Er{r}_{LSM}=\sum _{i=1}^{}{\left(x_i-y_i\right)}^2$$

(18)

where $x_i$ are the data points to be fitted, $y_i$ are the values of the corresponding solution and $Er{r}_{LSM}$ the regression error.

The fractional and conventional model parameters that resulted in the least regression error were calculated using an iteration algorithm [50,51] developed in Mathematica [35]. This same software was used to carry out the statistical analysis hereby described.

Additionally, a linear isotherm was considered to account for soil–solute interactions, resulting in the presence of a retardation factor R that divides the flow velocity $v^{*}=v/R$ and dispersion coefficients $D^{*}=D/R$ and ${\mathcal{D}}^{*}=\mathcal{D}/R$ [49,52]. This way, D stands for the dispersion coefficient in situations where no interaction takes place between the soil and the solute (ionic changes etc), and the diffusion process is solely governed by the soil’s porous network paths and the type of contaminant under consideration.

The independent variables x, t and v were obtained in the laboratory tests as the columns height (6 cm and 12 cm), test duration and average effluent flow velocity, respectively. For the average effluent flow velocity, a value of $4.9\times {10}^{-5}$ m/s and $10.58\times {10}^{-5}$ m/s for the samples with heights of 6 cm and 12 cm, respectively.

The retardation factor R for the experiments was calculated using the method described in Shackelford [53], i.e., the value in the $(V/V_v)$×$(C/C_0)$ breakthrough curve for which $C/C_0=0.5$, where V is the effluent volume and $V_v$ the voids volume.

Such an R estimation method is derived from the use of solution (5), which, given its infinite medium boundary condition, is symmetric around the point $C/C_0=0.5$ and allows for R to be calculated as the intercept of the $V/V_v$ breakthrough curve and the $C/C_0=0.5$ line. Therefore, for semi-infinite medium and thus asymmetrical solutions, i.e., Equations (3) and (17), a correction factor $f_c$ ($v^{*}=v{f}_c/R$) is used to further adjust the breakthrough curves and the flow velocity observed during the tests. The correction factor was not required for solution (6) due to its equivalency to solution (5) and its infinite medium boundary condition.

Figure 9 presents the fitting results for solutions (3), (6) and (17) with copper contamination and Figure 10 with lead contamination, both for the 6 cm soil column samples. In all the figures, the horizontal axis indicates the variable plotted (time) and its unit (days). Additionally, the code SH-XX has been used to name the samples, where H is the height of the sample and XX stands for the numbering sequence of the sample’s replica.

Table 2. Fitted parameters from 6 cm soil column samples with copper (${\mathrm{Cu}}^{2+}$) contamination.

	D	$\mathcal{D}$	$\mathit{\alpha }$	$\mathit{\beta }$	${\mathit{f}}_{\mathit{c}}$	R
	S6-01
Equation (3)	0.012	-	-	-	0.93	25.5
Equation (6)	-	0.684	1.22	-	-	25.5
Equation (17)	-	0.332	1.10	−0.88	$4.48\times {10}^{-9}$	25.5
	S6-02
Equation (3)	0.006	-	-	-	0.95	24.4
Equation (6)	-	0.056	1.52	-	-	24.4
Equation (17)	-	0.158	1.21	−0.86	0.58	24.4
	S6-03
Equation (3)	0.034	-	-	-	0.93	9.3
Equation (6)	-	0.420	1.46	-	-	9.3
Equation (17)	-	0.387	1.38	−0.79	0.83	9.2
	S6-04
Equation (3)	0.055	-	-	-	0.79	27.7
Equation (6)	-	1.38	1.30	-	-	27.7
Equation (17)	-	0.740	1.26	−0.98	$1.42\times {10}^{-9}$	27.7

and

Table 3. Fitted parameters from 6 cm soil column samples with lead (${\mathrm{Pb}}^{2+}$) contamination.

	D	$\mathcal{D}$	$\mathit{\alpha }$	$\mathit{\beta }$	${\mathit{f}}_{\mathit{c}}$	R
	S6-01
Equation (3)	0.009	-	-	-	0.95	36.6
Equation (6)	-	8.320	1.03	-	-	36.6
Equation (17)	-	0.328	1.11	−0.94	$6.84\times {10}^{-8}$	36.6
	S6-02
Equation (3)	0.008	-	-	-	0.93	35.5
Equation (6)	-	0.009	1.99	-	-	35.5
Equation (17)	-	0.276	1.15	−1.00	$2.08\times {10}^{-11}$	35.5

present the respective fitted transport parameters, where D is the dispersion parameter fitted by using Equation (3); $\mathcal{D}$ is the fractional dispersion parameter fitted by using Equations (6) and (17); $\alpha$ is the fractional derivative index; $\beta$ is the asymmetry parameter; $f_c$ is the correction factor (needed to adjust the breakthrough curves and flow velocity, as indicated) and R is the retardation factor (accounting for a linear isotherm of soil-contaminant interaction).

One can see from the regression errors $Er{r}_{LSM}$ in Figure 9 and Figure 10 that solution (17) presents the best fit to the experimental data and consequently generates the least regression error for each breakthrough curve. Additionally, the goodness-of-fit of solution (17) can be graphically observed in the results, especially in the tail portions of the breakthrough curves, where the asymmetry parameter $\beta$ allows an accurate simulation of the left and right tail portions of the curves. The fitted values of derivation index $\alpha$, which range closer to the lower limit of 1, give an indication that indeed an anomalous or super-dispersion of particles occurred in the column tests, especially considering they were conducted while preserving the original voids disposition of the soil in its natural state. The addition of the asymmetry parameter $\beta$ and its improvement in simulating the experimental data indicate that an anomalous and also preferential dispersion of particles took place, where fitted values of $\beta$ closer to the limit of $-1$ represent super-dispersion of particles in the opposite direction of flow. Such indicated behaviour may be a result of the upward flow of water used by the test equipment as opposed to the natural geological formation of the soil.

On the other hand, Figure 9c shows graphical proximity in the fitness of solution (6) and (17), although solution (17) resulted in a lower regression error. This may indicate that dispersion in soil sample S6-03 was anomalous but symmetrical. For soil sample S6-02 with lead contamination (Figure 10b), Equations (3) and (6) resulted in the same fit, while only the new solution in (17) was able to precisely fit the experimental data. The fit for S6-03 seems to indicate that a more complex interaction between the soil sample and the contaminant took place, which explains why the value of the retardation factor R was so different for this particular sample. In other words, the linear isotherm considered (which is represented by R) was not accurate for this sample. This is not the general case, as the other three replicated samples of S6 all behaved similarly.

The correction factor $f_c$ indicates how much the plume velocity $v^{*}=v{f}_c/R$ needed to be changed, due to the asymmetry of the model, for the minimization of the fitting error. Thus, it is observed that the conventional model (ADE) fits $f_c$ values close to one, indicating little difference in asymmetry considering the delay coefficient R value calculated for the symmetric solution. The new proposed fractional model, on the other hand, has asymmetry imposed by the boundary condition and the parameter $\beta$, resulting in $f_c$ values far from one, i.e., the plume velocity had to be modified to adjust the curves, except for sample S6-03, where the fit was similar to the conventional model. Finally, the fractional solution of Benson et al. [12] does not present a correction factor $f_c$ because it is symmetric and the adjustment of its plume velocity is done in the calculation of the delay factor R.

It is important to highlight that the transport experiment presented consisted of injecting a multi-species solution comprising nitrates and several heavy metals into contaminated soil. Therefore, the fate and transport of these chemical species is complex and may require more complicated reactive transport models than the one given by the FADE model with a simple retardation factor. This can be observed by noticing that not every contaminant BTC could be properly modelled.

Even though this apparent incompatibility shows up, the equations hereby developed were able to mathematically predict the overall behaviour of the system. In order to directly relate the parameters obtained to the chemical species studied, following works will extend the FADE model to multi-species reactive transport with a broader geochemical description of the processes involved.

Scale Effect

To understand the scale effect or scale-dependency of each model and solution, i.e., the precision to simulate particle dispersion at bigger scales with parameters obtained at smaller scales, the same column tests were conducted with 12 cm soil column samples.

Table 4. Fitted parameters from 12 cm soil column samples with copper (${\mathrm{Cu}}^{2+}$) contamination.

	D	$\mathcal{D}$	$\mathit{\alpha }$	$\mathit{\beta }$	${\mathit{f}}_{\mathit{c}}$	R	${\mathit{Err}}_{\mathit{LSM}}$
	S12-01
Equation (3)	0.043	-	-	-	0.95	11.9	0.038
Equation (6)	-	0.184	1.60	-	-	11.9	0.024
Equation (17)	-	0.136	1.61	−0.99	0.93	11.9	0.009
	S12-02
Equation (3)	0.066	-	-	-	0.95	8.1	0.021
Equation (6)	-	0.314	1.59	-	-	8.1	0.013
Equation (17)	-	0.234	1.59	−0.99	0.90	8.1	0.008
	S12-03
Equation (3)	0.029	-	-	-	0.97	13.3	0.072
Equation (6)	-	0.519	1.34	-	-	13.3	0.029
Equation (17)	-	0.235	1.39	−0.50	0.92	13.3	0.026
	S12-04
Equation (3)	0.992	-	-	-	0.53	23.6	0.039
Equation (6)	-	3.251	1.64	-	-	23.6	0.036
Equation (17)	-	8.916	1.18	−0.16	$4.92\times {10}^{-7}$	23.6	0.029

presents the fitted parameters for the 12 cm column samples with copper contamination and the respective LSM error. One can observe the same improvement in goodness-of-fit of the new solution (17) when compared to solutions (6) and (17), as shown by the reduction in regression error $Er{r}_{LSM}$.

In possession of the transport parameters fitted at both scales, one can statistically observe the change in the dispersion coefficient of each model. Figure 11 shows the statistical distribution of the relation between the dispersion coefficients fitted with 12 cm and 6 cm soil column samples, i.e., $D_{12cm}/D_{6cm}$ and ${\mathcal{D}}_{12cm}/{\mathcal{D}}_{6cm}$. The statistical analysis is composed of the 16 combinations between the dispersion coefficients of four 12 cm samples and four 6 cm samples, for each solution. A Student’s t-distribution was adopted for variance determination to account for the low data sample.

Table 5. Mean and variance of the increase in conventional and fractional dispersion coefficients when comparing fitted results from 12 cm and 6 cm soil column samples.

	${\mathit{D}}_{12\mathbf{cm}}/{\mathit{D}}_{6\mathbf{cm}}$		${\mathcal{D}}_{12\mathbf{cm}}/{\mathcal{D}}_{6\mathbf{cm}}$
	$\mathbf{\mu }$	${\mathbf{\sigma }}^{\mathbf{2}}$	$\mathbf{\mu }$	${\mathbf{\sigma }}^{\mathbf{2}}$
Equation (3)	21.74	2081.32	-	-
Equation (6)	-	-	5.94	198.41
Equation (17)	-	-	7.91	239.91

presents the mean $\mu$ and variance ${\sigma }^2$ values of the change in dispersion coefficient for each model.

From Table 5 and Figure 11, the conventional dispersion coefficient needs a mean increase of 22 times to simulate dispersion at a larger scale, while the fractional dispersion coefficients need a smaller increase of 6 to 8 times.

Such behaviour could indicate that the correct scaling effect of the fitted fractional diffusivity parameter may be physically justified by a fractal understanding of the porous media involved [54]. Assuming that the anomalous diffusive process scaling is modelled by the fractional differential equation index $\alpha$ and that the scaling effect of the fractal geometry of the medium can be modelled as a power law, then:

$$\mathcal{D}=l^{d_f}$$

(19)

where $\mathcal{D}$ represents the diffusivity coefficient of the medium, l is its characteristic length and $d_f$ is the fractal dimension of the medium.

For the present case, ${\mathcal{D}}_{12cm}/{\mathcal{D}}_{6cm}=2^{d_f}$. Since $2<d_f<3$ for soils in general, a ratio of $4<{\mathcal{D}}_{12cm}/{\mathcal{D}}_{6cm}<8$ is expected. This was experimentally obtained in the column tests and fitting analysis performed, which is a great insight that the fractional derivative solutions are adequate models for the phenomena observed.

5. Conclusions

The present work introduces a solution to FADE that satisfies a continuous source initial condition and semi-infinite boundary condition, which improves the evaluation of the initial concentration $C_0$ when close to the top of the soil layer, i.e., when $x=0$ and accounts for asymmetrical anomalous dispersion by considering the asymmetry parameter $\beta$. Furthermore, by using a Riesz fractional derivative definition, the R-FADE removes the $|cos(\pi \alpha /2)|$ factor from the solution, resulting in a scaling parameter $\gamma$ that does not quickly go to zero as $\alpha \to 1$.

Column tests at two different scales (6 cm and 12 cm) were conducted to verify the new solution and compare its performance to the conventional solution [4] and fractional solution given in Benson et al. [12]. The fitting of laboratory breakthrough curves with a non-linear LSM algorithm showed the improvement of the new solution to simulate the heavy-tailed curves, which is seen on the reduction in regression error $Er{r}_{LSM}$ and is graphically visible in Figure 9 and Figure 10.

Additionally, by statistically comparing the change in dispersion coefficients obtained at the smaller and larger scale, Figure 11 shows that the conventional model needs a wider increase in the value of the dispersion coefficient to simulate the larger scale. By considering a fractal power-law scaling function, the FADE and R-FADE solutions show their correctness of modelling not only the breakthrough curves, but also the fractal nature of the medium. Thus, the addition of a fractional index of derivation $\alpha$ and an asymmetry parameter $\beta$ to the contaminant transport model better describes the complexity of the dispersion of solutes in porous media and more accurately extrapolate laboratory tests to larger scales.