Numerical Simulation Of Nonlinear Advection-Diffusion-Fractional Equation of Pollutants in Porous Media

1. Introduction

The transport of dissolved contaminants through heterogeneous hydrogeology relies on the expression of their functional parameters. Non-linear uptake is a parameter that inherits significant ambiguities in its estimation when subjected to heterogeneous flow conditions [1]. The transport processes that govern the behavior of solutes in heterogeneous subsurface media are advection, dispersion, diffusion, sorption, and degradation. Among these processes, sorption and degradation play an important role in the infiltration and distribution of the solute in the heterogeneous porous medium [2]. The advection-dispersion equation (ADE) is the most commonly used method by researchers for the modeling of contaminant transport in subsurface environments. It describes contaminant transport processes in inhomogeneous porous media. However, it is less effective in describing this phenomenon in heterogeneous porous media [3,4,5]. EAD is based on the assumption of Fick's first law, whereas the solute transport process in heterogeneous porous media corresponds to a non-Fickian or anomalous mechanism [3]. Non-local transport models are alternative models proposed to explain the non-Fickian transport process. These non-local transport models include the spatial fractional advection-dispersion equation (s-FADE) and the temporal fractional advection-dispersion equation (t-FADE). The FADE equation is used to model the temporal and spatial variation of the concentration of one or more molecules of a substance. The fractional advection-dispersion-reaction model plays a key role in studying non-linear physical problems [6]. In groundwater hydrology, fractional advection-dispersion equations represent the passive transfer of tracers in the flow of a porous fluid [7] . Many physical models used in engineering and science are modeled using fractional derivatives. Over the last few decades, the use of fractional derivatives by researchers has increased dramatically in many important fields, including chemistry, signal theory, information systems, economics, viscoelasticity, finance, physics, and fluid mechanics [8]. Applications of fractional calculus have been successful in modeling various types of fractional differential equations, both linear and non-linear, Prashant Pandey et al [9,10]. Fractional derivatives are applied to many physical problems associated with dynamical systems involving differential equations [11]. Transport differential equations based on fractional derivatives are particularly well suited to describing anomalous dispersion processes, where a plume of particles propagates at a velocity incompatible with the classical Brownian motion model. When a fractional derivative replaces the second derivative in a diffusion or dispersion model, diffusion is improved [12]. In the groundwater flow model, the second spatial derivative in the traditional FADE is replaced by a fractional derivative of order α, where 1 < α < 2 [13]. FADEs have been successfully applied to simulate solute transport in homogeneous and heterogeneous, saturated and unsaturated soils on a laboratory scale [14], in a sand and gravel aquifer [5,15] and in a highly heterogeneous aquifer [15]. The approximate numerical solution of the reaction-diffusion-fraction equation remains a focus of important research today because of its nature and application to petroleum reservoir simulation modeling, fluid flow in porous media, global water production, power generation, and drinking water production [6] .

The use of fractional advection-dispersion models has been discussed by many researchers from different angles. Here are some examples. [8] used the energy method to show the unconditional stability and order of convergence obtained from the discrete-time scheme of the s-FADE in the finite domain, where the time and space derivatives are the fractional Caputo derivative. [16] studied the existence and uniqueness of the fractional advection-dispersion equation described by the generalised Caputo-Liouville fractional derivative. The method of Laplace transformations was developed to obtain approximate solutions of the fractional advection-dispersion equation. In [17], the combination of an advection-dispersion equation in fractional time and fractal space was studied using fractional fractal derivatives to model groundwater transport. In the same work, [17] showed that the fractal advection-dispersion equation shows that the fractional order in time controls the peak location of the breakthrough curve and the fractal dimension in space controls the drag in the simulated breakthrough curve. The researcher [18] developed a numerical and experimental method to solve the fractional advection-dispersion equation as a relation used in groundwater by exploring the applicability developed by the numerical method of homogeneous and heterogeneous solutions in subsurface media. The authors of the study [19] developed a fractional dispersion aversion model through experimental and numerical studies to evaluate the solute transport parameter of s-FADE in homogeneous and heterogeneous soils.

Non-linear fractional partial differential equations are not as easy to solve as linear equations. Furthermore, there is no exact solution to the fractional order differential equation by any method. Therefore, the fractional-order advection-dispersion fractional equation (FADE) makes research more important and practical [7] . It can more accurately describe the transport of contaminants in unsaturated porous media.

A mathematical model with four types of input conditions was developed by [1] to evaluate the variability of the dispersion coefficient at temporal and spatial scales on the contaminant transport mechanism in unsaturated porous media subjected to a non-linear sorption model. The variability of the sorption parameters of the Freundlich model on contaminant dispersion in aquifers was not considered by these authors.

In our recent work [20,21,22], the numerical solution of the classical dispersion and advection equation has been developed taking into account the time and space dependent dispersion coefficient to describe the contaminant transport mechanism in subsurface media. The latter does not consider the fractional advection-dispersion equation (FADE) model, which is now very effective in adequately capturing contaminants in heterogeneous porous media [4].

[23] Rashmi Radha et al (2023) numerically solved a two-dimensional contaminant transport model with time-varying axial input sources subject to nonlinear sorption, decay and production. They highlighted the influence of hydrogeological parameters such as medium porosity, density, sorption conditions and dispersion coefficients, but failed to highlight the influence of the adsorption coefficient (Kd) on the spatio-temporal variation of contaminant concentration in subsurface media. Furthermore, these studies did not consider the fractional advection-dispersion equation (FADE) model, which remains essential for describing the level of contamination in non-linear porous media.

This paper proposes a numerical solution of the fractional advection-dispersion equation (FADE) model to explicitly evaluate and simulate the contaminant transport mechanism in unsaturated porous media subjected to a non-linear sorption model.

We aim to evaluate the influence of parameters n (sorption intensity or Freundlich parameter) and adsorption coefficient (Kd) on the spatiotemporal variation of contaminant concentration in the underground environment. This work will also make it possible to select the values of n (sorption intensity) and the adsorption coefficient (Kd) at the aquifer outlet to optimize the retention of pollutants in the underground environment.

2. Materials and Methods

2.1. Physical Model of System

Water containing pollutants, such as sewage and radioactive waste, seeps into the soil matrix as does direct runoff from surfaces resulting from runoff. This water eventually enters the aquifer and is influenced by distance and subsoil parameters, and becomes a source of drinking water. As water passes through the soil, pollutants are mixed, adsorbed, dispersed, and diffused in the flow. The classic porous domain considered for the present study is shown in Figure 1 below.

2.2. Mathematical Reformulation

Pollutants infiltrating groundwater can be studied by integrating zero-order and first-order adsorption and decay coefficients. The movement of solute in the flow region occurs due to advection and dispersion. One-dimensional solute motion in a heterogeneous porous medium can be studied numerically using the spatial fractional advection-dispersion equation for a one-dimensional unsteady-state velocity field with a first-order decay function, and an absorption coefficient non-linear. s-FADE is a spatially nonlocal model capable of describing solute particle motions encompassing a very large transition due to high heterogeneity to preferential high-velocity flow paths. The influence of heterogeneity leads to the variability of some coefficients of the FADE model hence the fractional-order α in this work is allowed to be variable in the porous medium, a fraction instead of the integer constant of 2. The one-dimensional advection-dispersion space-fractional equation for a conservative contaminant can be written as follows [1,18]:

$$R_d\frac{\partial c(x,t)}{\partial t}=D_f\frac{{\partial }^{\alpha }c(x,t)}{\partial x^{\alpha }}-v_{0x}\frac{\partial c(x,t)}{\partial x}-\lambda c(x,t)$$

(1)

Where c is the concentration of the contaminant, $D_f$ is the fractional dispersion coefficient, t is the time, and x is the distances traveled by the fluid in the porous medium. λ is the first-order decay coefficient and $R_d$ is the delay factor generated by the nonlinear sorption of the pollutant. The temporal fractional advection-diffusion model can be written in the form of the following equation (3) [1,17]:

$$\frac{{\partial }^{\gamma }c(x,t)}{\partial t^{\gamma }}=\frac{1}{R_d}\left[D_f\frac{\partial c(x,t)}{\partial x}-v_{0x}\frac{\partial c(x,t)}{\partial x}-\lambda c(x,t)\right]$$

(2)

2.3. Description of the Coefficient of the First Order Decay Coefficient

A first-order equation is used to represent the kinetics of simple reactions such as certain microbial degradation reactions. The equation 3 is written as:

$$\frac{\partial c}{\partial t}=-\lambda c$$

(3)

Where c is the concentration in the liquid phase, λ is the degradation coefficient which is expressed as a function of the half-life of the compound:

$$\lambda =\frac{\ln (2)}{T_{1/2}}$$

(4)

Where $T_{1/2}$ is the half-life time. The degradation kinetics of organic pollutants strongly depend on the properties of the medium (presence of an electron sensor, such as oxygen). Moreover, the kinetics observed for certain conditions is not necessarily transposable to other conditions, [24].

2.4. Description of the Model of the Delay Factor Generated by the Sorption of the Pollutant (Rd)

There are two types of sorption isotherm: linear and non-linear. The linear sorption isotherm applies to low concentrations of contaminants and to absorption. If the sorption isotherm is linear, there is no upper limit to contaminant sorption to soil particles. However, solute transport is more realistic for the non-linear sorption model [23] . The Freundlich transport model has shown that below a threshold concentration, the retardation factor increases with the Freundlich isotherm exponent n, whereas the trend is reversed and the retardation factor decreases as n increases above this specific concentration. Due to hydrodynamic dispersion, the concentrations in the plume decrease towards the edges of the plume. The model for the non-linear absorption dependent retardation factor is given by the relationship:

$$R_d=1+\frac{{\rho }_b}{\theta }{k}_{d}n{c}^{n-1}$$

(5)

Kd is the adsorption coefficient and n is the sorption intensity. (${\mathit{\rho }}_{\mathit{b}}$) is the specific weight of the solid particle in the matrix, $\theta$ is the porosity of the porous medium (The porosity is the ratio between the volume of all the pores of a sample of the soil representative of the medium considered, and the total volume of the same sample). $k_d$ expresses the ratio between the amount of solute concentration absorbed per unit area and the amount of concentration in solution. It eventually causes a delay of solute relative to mobile fluid in the fracture. The delay factor increases with the Freundlich isothermal exponent n. The retardation coefficient increases non-linearly with decreasing concentration. The effect of the isothermal exponent on the lag coefficient depends on the order of magnitude of the concentration [9].

2.5. Initial and Boundary Conditions

From the source of the pollutants, the concentration disperses further into the porous medium in the direction of flow (x) Figure 1. A one-way horizontal flow of water through an aquifer of finite length L is considered [21]. The problem can be mathematically defined by equations 1 and 2 with the initial and boundary conditions given successively in equations (3) and (4) as follows:

$$c(x,0)=0,0\prec x\le L,$$

(6)

$$\left\{\begin{array}{l}c(0,t)=c_0,\forall t\ge 0\\ \frac{\partial c(L,t)}{\partial x}=0,\forall t\succ 0\end{array}\right.,$$

(7)

2.6. Numerical Solution of the Mathematical Model

The one-dimensional solute transport equation is a parabolic type partial differential equation in which, the finite difference technique is commonly used to obtain the numerical solution. An explicit finite difference technique is generally used and eventually gives rise to the restrictive stability criteria [10,20].

2.7. Spatial Fractional-Order Derivative

In this part, the basic theory of the numerical solution of the fractional space of equation 1 in an infinite domain is [18]. The discrete approximation of the fractional derivative term can be defined from the Grunwald norm, as follows.

$$\frac{{\partial }^{\alpha }c}{\partial x^{\alpha }}=\frac{1}{{\mathsf{\Gamma }}_{(-\alpha )}}\underset{n\to \infty }{\lim }\frac{1}{h^{\alpha }}{\displaystyle \sum _{k=0}^i\frac{{\mathsf{\Gamma }}_{(k-\alpha )}}{{\mathsf{\Gamma }}_{(k+1)}}}c(x-kh,t),$$

(8)

Where i is a positive index, ${\mathsf{\Gamma }}_{(-\alpha )}$ is the gamma function. The value of the fractional derivative at a point x depends on the values ​​of the function at that point and all points to the left of that point of interest x. The explicit Euler method associated with the Grunwald approximation on fractional derivatives will be used in this part, which is as follows:

$$\frac{c_i^{j+1}-c_i^j}{\mathsf{\Delta }t}=c_i^j+\frac{1}{{(R_d)}_i^j}\left[-v_{0x}\frac{c_i^j-c_{i-1}^j}{h}+\frac{D_f}{h^{\alpha }{\mathsf{\Gamma }}_{(-\alpha )}}{\displaystyle \sum _{k=0}^i\frac{{\mathsf{\Gamma }}_{(k-\alpha )}}{{\mathsf{\Gamma }}_{(k+1)}}}{c}_{i-k}^j\right]$$

(9)

were (i= 1,2,..., $N_x$), (k=0, 1,2, ...,) et (j=0, 1,2, ..., $N_t$) $x=ih,t=j\mathsf{\Delta }t$ . Normalized Grunwald weights are defined by:

$$f_k=\frac{{\mathsf{\Gamma }}_{(k-\alpha )}}{{\mathsf{\Gamma }}_{(-\alpha )}{\mathsf{\Gamma }}_{(k+1)}}$$

(10)

These normalized weights only depend on the order α and the index k. (For example, the first four terms of this sequence are given by :$f_0$=1, $f_1$=$-\alpha$, $f_2$= $\alpha (\alpha$-1)/2!, $f_3$= $-\alpha (\alpha$-1) $(\alpha$-2)/3!, the resulting equation to be solved explicitly is reduced by:

$$c_i^{j+1}=\frac{1}{{(R_d)}_i^j}\left[\left({(R_d)}_i^j+\lambda -\frac{\mathsf{\Delta }t}{h}{v}_{0x}+\frac{\mathsf{\Delta }t}{h^{\alpha }}{D}_f\right)c_i^j+\left(\frac{v_{0x}}{h}+\frac{\alpha }{h^{\alpha }}{D}_f\right)\mathsf{\Delta }t{c}_{i-1}^j+\frac{D_f}{h^{\alpha }}{\displaystyle \sum _{k=2}^i{f}_k}{c}_{i-k}^j\right]$$

(11)

2.8. Temporal Fractional-Order Derivative

In this section, the solution of the temporal fractional numerical equation (2) in an infinite domain is developed. The Caputo-Liousville fractional derivative is exploited with the following equation Lui et al, (2012):

$$\frac{{\partial }^{\gamma }c}{\partial t^{\gamma }}=\frac{1}{\mathsf{\Delta }{t}^{\gamma }{\mathsf{\Gamma }}_{(2-\gamma )}}{\displaystyle \sum _{k=0}^j{w}_k^{\gamma }}\left[c(x,t_{j+1-k})-c(x,t_{j-k})\right],$$

(12)

were $w_k^{\gamma }={(k+1)}^{1-\gamma }-k^{1-\gamma },k=1,2,\mathrm{...},n$ ,(i=0, 1,2,..., Nx) et (j= 1,2, ..., N_t) $x=ih,t=j\mathsf{\Delta }t$ . The explicit Euler method associated with the Caputo-Liousville approximation on the fractional derivatives of equation (2) is used in this part, which is presented by the following equation (13):

$$\frac{{(R_d)}_i^j}{\mathsf{\Delta }{t}^{\gamma }{\mathsf{\Gamma }}_{(2-\gamma )}}{\displaystyle \sum _{k=0}^j{w}_k^{\gamma }}\left[c(x,t_{j+1-k})-c(x,t_{j-k})\right]=-\frac{v_{0x}}{h}\left(c_i^{j+1}-c_{i-1}^j\right)+\frac{D_f}{h^2}\left(c_{i+1}^j-2c_i^j+c_{i-1}^j\right)+\lambda c_i^j$$

(13)

By simplifying equation (13) we obtain the following equation (14) to be solved explicitly:

$$c_i^{j+1}=c_i^j-a_1{\displaystyle \sum _{k=1}^j{w}_k^{\gamma }}\left[C(x,t_{j+1-k})-c(x,t_{j-k})\right]-a_2\left(c_i^j-c_{i-1}^j\right)+a_3\left(c_{i+1}^j-2c_i^j+c_{i-1}^j\right)$$

(14)

Avec $a_1=w_0^{-\gamma }$,$a_2=\frac{w_0^{-\gamma }{v}_{0x}}{b\mathsf{\Delta }t}$, $b=\frac{{(R_d)}_i^j\mathsf{\Delta }{t}^{-\gamma }}{{\mathsf{\Gamma }}_{(2-\gamma )}}$,$a_3=\frac{D_i{w}_0^{\gamma }{h}^{-2}}{b}$,${(R_d)}_i^j=1+\frac{{\rho }_b}{{\theta }_m}{k}_{d}n{(c^{n-1})}_i^j$

The discretization of the initial conditions and limit associated with the equation (6) and (7) to apply to this method is given successively to the following equation (15) and (16):

$$c_i^0=0,0\prec i\le k-1,$$

(15)

$$\left\{\begin{array}{l}{c}_0^j=c_0,\forall j\ge 0\\ \frac{\partial c_{N_x}^j}{\partial x}=0,\forall j\ge 0\end{array}\right.,$$

(16)

2.9. Stability Condition for the Spatial Fractional Order Derivative

Suppose that $c_i^0$ is the only term that has an error, so the perturbed value is $\underset{¯}{c_i^0}=c_i^0+{\delta }_i^0$ , this perturbation produces a perturbed value for given $\underset{¯}{c_i^1}=c_i^1+{\delta }_i^1$ , equation (11) becomes:

$$\underset{¯}{c_i^1}={\beta }_i{c}_i^0+\frac{1}{{(R_d)}_i^0}\left[\left(\frac{v_{0x}}{h}+\frac{\alpha }{h^{\alpha }}{D}_f\right)\mathsf{\Delta }t{c}_{i-1}^j+\frac{D_f}{h^{\alpha }}{\displaystyle \sum _{k=2}^i{f}_k}{c}_{i-k}^j\right]={\beta }_i{\delta }_i^0+c_i^1$$

(17)

where the factor

$${\beta }_i=\frac{1}{{(R_d)}_i^0}\left[{(R_d)}_i^0+\lambda -\frac{\mathsf{\Delta }t}{h}{v}_{0x}+\frac{\alpha }{h^{\alpha }}{D}_f\right]$$

(18)

So we have ${\delta }_i^1={\beta }_i{\delta }_i^0$. That is, the error is magnified by the factor ${\beta }_i$ when the finite difference equation is advanced one time step. After j no time, we can write ${\delta }_i^j={\beta }_i^j{\delta }_i^0$ which is the amplification factor. For the explicit Euler method to be stable, $\left|{\beta }_i\right| \le 1$ or

$$\frac{\mathsf{\Delta }t{v}_{0x}}{h}-(\lambda +1+R_d{}_i)\le \frac{\alpha D_f}{h^{\alpha }}\le \frac{\mathsf{\Delta }t{v}_{0x}}{h}-(\lambda +1-R_d{}_i)$$

(19)

2.10. Stability Condition for the Temporal Fractional Order Derivative

In this part the equation (14) can still be written in the form:

$$c_{i+1}^j=\left(2+\lambda -\frac{1}{a_3}+\frac{a_2}{a_3}\right)c_i^j-\left(1+\frac{a_2}{a_3}\right)c_{i-1}^j+\frac{1}{a_3}{c}_i^{j+1}+\frac{a_1}{a_3}{\displaystyle \sum _{k=1}^j{w}_k^{\gamma }}\left[c(x,t_{j+1-k})-c(x,t_{j-k})\right]$$

(20)

Suppose that $C_0^j$ is the only term that contains the error, so the perturbed value is $\underset{¯}{c_0^j}=c_0^j+{\delta }_0^j$This perturbation gives rise to another perturbed value $C_1^j$which is $\underset{¯}{c_1^j}=c_1^j+{\delta }_1^j$, then Equation (20) becomes:

$$\underset{¯}{c_1^j}=m_i{c}_0^j+\frac{1}{a_3}{c}_i^{j+1}-\left(1+\frac{a_2}{a_3}\right)c_{i-1}^j+\frac{a_1}{a_3}{\displaystyle \sum _{k=1}^j{w}_k^{\gamma }}\left(c_i^{j+1-k}-c_i^{j-k}\right)=m_i{\delta }_0^j+c_0^j$$

(21)

Where the factor $m_i=2+\lambda -\frac{1}{a_4}+\frac{a_2}{a_4}$ is the amplification factor associated with the error with ${\delta }_1^j=m_i{\delta }_0^j$ when the finite difference of the equation is advanced by a spatial step, for i spatial step , we have ${\delta }_i^j=m_i^j{\delta }_0^j$ . For Euler's method to explain is stable, $\left|m_i\right| \le 1$ which boils down to the following relationship (22):

$${\mathsf{\Gamma }}_{(2-\gamma )}\left(\frac{D_f\mathsf{\Delta }t-v_{0x}{h}^2}{{(R_d)}_0^j{w}_0^{\gamma }}\right)\le \frac{h^2}{\mathsf{\Delta }{t}^{1+\gamma }}\le {\mathsf{\Gamma }}_{(2-\gamma )}\left(\frac{3D_f\mathsf{\Delta }t-v_{0x}{h}^2}{{(R_d)}_0^j{w}_0^{\gamma }}\right)$$

(22)

The experimental values of the delay factor parameters generated by the non-linear sorption have been exploited in the work of [2], the same data adding to the data used in the work of [1] will be exploited in this work, see Table 1 below

Table 1. Parameter data of delay factor R parameters.

Table 1. Parameter data of delay factor R parameters.

3. Results:

3.1. Effect of Spatial Fractional Order on the Concentration of Pollutants in a Porous Medium

In this section, the discrete numerical solution of the equation (11) associated with the initial conditions and the discrete limits of the equations (15) and (16) is simulated to obtain a certain number of results.

Euler's explicit method associated with Grunwald's fractional derivative approximation was used to simulate the numerical solution of equation (8) as a function of aquifer length, time, constant dispersion coefficient, and the flow rate of fluids in the aquifer. For this purpose, the parameters of the model used were selected in previous studies [1,2,14,19,20]. Figure 3 presents the Spatial variation of the concentration of pollutants in the porous medium for different spatial fractional-order values ​​α (α =1.8, 1.85, 1.9, 1.95.2) for the absorption coefficient Kd=0.75 L/g and n=0.7. Figure 3 shows that the mass of solute retained in the underground medium increases with (α). The concentration profiles follow an exponential decrease in the aquifer independently of the spatial fractional order (α). This result is similar to those of [1] for the case of the classic ADE model where α=2. These authors numerically simulated a mathematical model with four types of input conditions, to assess the variability of the dispersion coefficient on temporal and spatial scales on the contaminant transport mechanism ants in saturated porous media. They did not take into account the spatial fractional advection-dispersion equation (s-FADE) model, which remains today the most complex model to assess the level of a contaminant in underground environments, which is the object of our work. It is also observed that the slope of the concentration profile becomes closer and closer with a decrease in the value of α.

Figure 2. Geometry of the problem.

Figure 2. Geometry of the problem.

Four points were selected in the aquifer (L/4, L/2, 3L/4, 0.9L) to assess the dispersion behavior of pollutants in the system. Figure 4 and Figure 5 show the relationship between the evolution of the concentration of pollutants at different points as a function of the absorption coefficients (Kd) and the absorption intensity (n) by successively using the equation model of fractional advection-dispersion (FADE) and the classical advection-dispersion equation (ADE) model. The highest values ​​of pollutant concentrations are obtained near the source (x = L / 4) while the low values ​​were found near the outlet of the system (x = 0.9 L). These results show that the dispersion of pollutants increases with the transport distance. The results also suggest that absorption coefficients (Kd) have a strong effect on pollutant dispersion rather than absorption intensity (n). We observe in Figure 5 a better retention of the concentration of pollutants due to the different degrees of heterogeneity generated by the model (FADE) to explain the process of solute transport in natural porous media. Figure 4 and Figure 5 give us information on the zones that were exploited by [2] from the entrance to the exit of the aquifer because in their work they took the value Kd=0.345 L/g and n=0.32. [1] meanwhile, he rather set the value Kd=0.75 L/g and n=0.7 which remains localized according to our study in the areas of average concentrations. In the light of Figure 4 and Figure 5, it emerges that the data of the absorption coefficients Kd and the intensity of absorption n exploited by their work are well located in a precise zone in the aquifer according to the selected points. For proper environmental modeling of solute transport in porous media, experimental studies should be performed. By specifying which point of the aquifer is exploited in relation to the source of the pollutant, to find the values ​​of absorption Kd and the intensity of absorption n, rather than fixing them as is done in the Natarajan literature [1]. It is also possible to find the range of these parameters that could improve system performance.

Figure 6 and Figure 7 show the evolution of pollutant concentrations as a function of different absorption constants (Kd) and absorption intensity (n) near the aquifer outlet (0.9L) for the two models s-FADE and ADE. These two figures allow us to find the range of these parameters Kd and n which can improve the performance of the system. It has been noted in Figure 6 that, for Kd equal to 0.2 L/g and all the values of n ∈R + in the aquifer, the lowest concentration values are observed with a range between 0 and 3% of the initial concentration, it has also been noted in Figure6 that, for Kd equal to 1 L/g and all the values of n ∈R +, the highest concentration values are observed with a range between 0 and 0.68 of the initial concentration. When analyzing the absorption constants (Kd) and the absorption intensity (n) for the model concentration profiles FADE, it appears that the appropriate values of Kd are less than or equal to 0.2L/ g for all values of n ∈R+. The high concentration values are obtained for Kd greater than or equal to 1L/g. Figure 6 shows that, for (Kd) equal to 1 L/g, the same trend is observed as in Figure 6 for the FADE model. However, for Kd equal to 0.2 L/g and all the values of n ∈R + in figure 7 for the ADE model, the lowest concentration values are observed with a range between 0 and 10% Co rather than 0 and 3% Co as observed in Figure 6. These values of Kd found are used at different points of the aquifer in this work to analyze the effect of the nonlinear absorption coefficient on the transport of contaminants in porous media.

Figure 8a,b successively illustrate the spatial distribution of the concentration of pollutants with selected constants of the absorption coefficient Kd = 0.2 L/g and Kd=0.5L/g. The concentration profiles according to space decrease exponentially in the aquifer for the values ​​of n (0.1, 0.4, 0.7) for the two models s-FAED and AED whatever the value of kd. This profile decreases rapidly for the model FAED, reflecting greater retention of the mass of solute in the medium. These results are similar and more complex to those obtained by [14] because they did not take into account the influence of absorption coefficients (Kd) and absorption intensity (n) in their work. These authors have even an experimental study associated with the model s-FADE to model the conservative solutes, in the end, their results indicate that the best simulation was obtained from the s-FADE model compared to the ADE model. The different concentrations all decrease together from the inlet of the aquifer independently of the model, to reach zero horizontal asymptotes for low absorption coefficients (Kd) figure 8a. And horizontal asymptotes at the different positions X = 0.29, 0.6, 0.8 respectively for three values ​​of n (0.2, 0.5, 0.7) for the s-FADE model when kd=0.5l/g figure 8a. The residual concentration percentage associated with these positions is X = 0.05; 0.08; 0.1. On the other hand, the different horizontal asymptotes are obtained rather at the different positions X = 0.5, 0.65, 0.9 for the ADE model respectively for the same values ​​of n (0.2, 0.5, 0.7). The associated residual concentration percentage is 0.1; 0.14; 0.16. These positions increase with the increase of n. The residual solute concentration percentage remains low for the s-FADE model. These results confirm that the increase in n assigned to the s-FADE models significantly reduces the mass of solute in the underground medium compared to the ADE models. This is hard to the different degrees of heterogeneity engendered by the model FADE to explain the process of solute transport in natural porous media [19].

3.2. Effect of Temporal Fractional Order on the Concentration of Pollutants in a Porous Medium

In this section, the discrete numerical solution of the equation (14) associated with the initial conditions and the discrete limits of the equations (15) and (16) is simulated to obtain the results.

Figure 9 shows the influence of temporal fractional order (γ) on the distribution of solute concentrations in a subterranean medium. The values ​​of the temporal fractional order are chosen in the order γ (γ = 0.2, 0.4, 0.6, 0.8, 1). Figure 9 further shows the impact of the absorption coefficient Kd (Kd = 0.2 and 0.75 l/g) on ​​the dispersion of pollutants in the environment. We note that there is low retention of the mass of solute when the value of the time-fractional order (γ) increases in the porous medium and approaches the source of the pollutants when the value of the coefficient of Kd absorption is low in the same medium, which contradicts the work of [26,27]. These authors have developed through the method of boundary elements the solution of the two-dimensional anomalous diffusion equation, and show the ability and necessity of small-time steps to capture the large temporal gradients that arise in the solution of problems governed by ADEs containing the fractional derivative. In their work, they did not allude to the parameter of the absorption coefficient Kd which plays the role of possibly causing a solute delay with respect to the mobile fluid in the underground medium. By analyzing these results, it appears that the retention of the mass of solute increases with absorption coefficient Kd. The advection dispersion model becomes classic when γ=1.

3.3. Validation Model

An Experimental solution on the contaminant transport equation in an underground environment with a constant source of contaminant provided by [18] is used as shown in the Figure 10 below. It is observed from Figure 10 that there is a good agreement between the numerical solution and the experimental solution from the model presented. The latter highlights the use of our model to describe and understand the level of contaminants in porous media.

4. Conclusion

In this study, the effect of the nonlinear absorption-dependent delay factor on solute transport in porous media was analyzed using two dispersion models ADE and FADE. The explicit Euler method associated with the Grunwald and Caputo-Liousville approximation on fractional derivatives was used to model the solute transport. When analyzing the results obtained, the following conclusions were drawn from this study:

• The amplitude of the initial concentration of pollutants is influenced by the different points taken in the aquifer for two dispersion models ADE and FADE used in the aquifer, and thus, the present models project that the resulting mobility of the contaminant from the groundwater is very sensitive to the boundary of the aquifer with which it is associated.
• A better solute mass retention is observed when the model is FADE due to the different degrees of heterogeneity generated by this model to explain the process of solute transport in natural porous media
• The sorption partition coefficient (Kd) plays a very significant role in the dispersion of contaminants in underground environments, regardless of the model used.
• The behavior of the concentration profiles remains the same for the different magnitudes of the absorption intensity (n) in the whole aquifer for the two dispersion models.