On the Generalized Mass Transport Equation to the Concept of Variable Fractional Derivative

Abstract

The hydrodynamic dispersion equation was generalized using the concept of variational order derivative. The modified equation was numerically solved via the Crank-Nicholson scheme. The stability and convergence of the scheme in this case were presented. The numerical simulations showed that, the modified equation is more reliable in predicting the movement of pollution in the deformable aquifers, than the constant fractional and integer derivatives.

1. Introduction

Anomalous dispersion phenomena are systematically pragmatic in physics, chemistry, and biology fields [1–4]. To distinguish irregular diffusion phenomena, constant-order fractional diffusion equations were initiated and have received breathtaking success. Nonetheless, it has been recognized that the constant order fractional dispersion equations are not competent of typifying some multifaceted dispersion processes, for example, dispersion process in inhomogeneous or heterogeneous medium [5]. In addition, when we think about dispersion process in porous medium, if the medium structure or external field changes with time, in these circumstances, the constant-order fractional diffusion equation model cannot be used to well illustrate such phenomenon [6, 7]. This is the case of the groundwater pollution problem; the medium through which the pollution is dispersing is heterogeneous and change with time. Motionless in some environmental science dispersion processes, the concentration of particles will determine the dispersion pattern [8, 9]. To solve the above problems, the variable-order (VO) fractional dispersion equation models have been suggested for use [10].

The leading edge work of VO operator can be traced to Samko et al. by introducing the variable order integration and Riemann-Liouville derivative in [10]. It has been acknowledged as a prevailing modelling approach in the fields of viscoelasticity [10], viscoelastic deformation [11], viscous fluid [12], and anomalous transmission [13]. The problem encounter in groundwater is that the geometry of the aquifer in which the flow or the pollution takes place is not well known. In addition, the geological formation through which the contamination takes place changes in time and space. These phenomena cannot therefore accurately be described neither via the classical hydrodynamic dispersion equation nor via the fractional version of hydrodynamic dispersion equation. In this paper, we generalize the hydrodynamic dispersion equation using the concept of the variable order derivative.

2. Definition and Problem Formulation

The fractional operators (fractional derivatives and integrals) refer to the differential and integral operators of arbitrary order, and fractional differential equations refer to those containing fractional derivatives. The former are the generalization of integer-order differential and integral operators and the latter are the generalization of differential equations of integer order. The operators of variable-order, which fall into a more complex category, are the derivatives and integrals whose orders are the functions of certain variables.

2.1. Variational Order Operator

Definition 1. Left and right Riemann-Liouville integrals of variable order: let 0<𝛼(𝑥,𝑡)<1 for all (𝑥,𝑡)∈[𝑎,𝑏] and 𝑓∈𝐿₁[𝑎,𝑏]; then, 𝑎𝐼_𝑡^{𝛼(⋅,⋅)}𝑓(𝑡)=^𝑡_𝑎1Γ[𝛼(()𝑥,𝑡)]𝑥−𝑡^{𝛼(𝑥,𝑡)−1}𝑓(𝑥)𝑑𝑥,(𝑡>𝑎) (1) is called the left Riemann-Liouville integral of variable fractional order 𝛼(⋅,⋅), while 𝑡𝐼_𝑏^{𝛼(⋅,⋅)}𝑓(𝑡)=^𝑏_𝑡1Γ[𝛼(()𝑥,𝑡)]𝑡−𝑥^{𝛼(𝑥,𝑡)−1}𝑓(𝑥)𝑑𝑥,(𝑡<𝑏) (2) is referred to as the right Riemann-Liouville integral of variable fractional order 𝛼(⋅,⋅).

Definition 2. Left and right Riemann-Liouville derivatives of variable fractional order: let 0<𝛼(𝑥,𝑡)<1 for all (𝑥,𝑡)∈[𝑎,𝑏]. If 𝑎𝐼_𝑡^{1−𝛼(⋅,⋅)}𝑓∈𝐴𝐶[𝑎,𝑏], then, the left Riemann-Liouville derivative of variable fractional order 𝛼(⋅,⋅) which is given as 𝑎𝐷_𝑡^{𝛼(⋅,⋅)}𝑓(𝑡)=𝑑𝑑𝑡^𝑡_𝑎1Γ[(()1−𝛼𝑥,𝑡)]𝑡−𝑥^{−𝛼(𝑥,𝑡)}𝑓(𝑥)()𝑑𝑥,𝑡>𝑎 (3) is called the left Riemann-Liouville derivative of variable fractional order 𝛼(⋅,⋅), while in the same line of idea we have the following expression: 𝑡𝐷_𝑏^{𝛼(⋅,⋅)}𝑓(𝑡)=𝑑𝑑𝑡^𝑏_𝑡1Γ[(()1−𝛼𝑥,𝑡)]𝑥−𝑡^{−𝛼(𝑥,𝑡)}𝑓(𝑥)()𝑑𝑥,𝑡<𝑏 (4) which is referred to as the right Riemann-Liouville derivative of variable fractional order 𝛼(⋅,⋅).

Definition 3. Left and right Caputo derivatives of variable fractional order: let 0<𝛼(𝑥,𝑡)<1 for all (𝑥,𝑡)∈[𝑎,𝑏]. If 𝑎𝐼_𝑡^{1−𝛼(⋅,⋅)}𝑓∈𝐴𝐶[𝑎,𝑏], then, the left Riemann-Liouville derivative of variable fractional order 𝛼(⋅,⋅) is given as 𝑎𝐷_𝑡^{𝛼(⋅,⋅)}𝑓(𝑡)=^𝑡_𝑎1Γ[(()1−𝛼𝑥,𝑡)]𝑡−𝑥^{−𝛼(𝑥,𝑡)}𝑑𝑓(𝑥)()𝑑𝑥𝑑𝑥,𝑡>𝑎 (5) which is called the left Caputo derivative of variable fractional order 𝛼(⋅,⋅), while in the same line of idea we have the following expression: 𝑡𝐷_𝑏^{𝛼(⋅,⋅)}𝑓(𝑡)=^𝑏_𝑡1Γ[(()1−𝛼𝑥,𝑡)]𝑥−𝑡^{−𝛼(𝑥,𝑡)}𝑑𝑓(𝑥)()𝑑𝑥𝑑𝑥,𝑡<𝑏 (6) which is referred to as the right Caputo derivative of variable fractional order 𝛼(⋅,⋅).
However, this work we will throughout use the following definition.

2.2. Problem Formulation

The concern here is the generalization of advection dispersion equation by including a possible effect of heterogeneity or variability of the aquifer into the mathematical formulation. A one-dimensional model consisting of an infinitely long homogeneous isotropic porous medium with a steady state uniform flow with a seepage velocity 𝑣 is considered here. A particular chemical is injected from one end of the model for a period of time 𝑡₀ such that the input concentration varies as an exponential function of time. The value of that chemical concentration at any time 𝑡 and at a distance 𝑥 from the injection boundary, allowing for the decay and adsorption, may be obtained from the solution of the following set of equations [14]; more details for this model can be found in [15–18]. Consider 𝐷𝜕²𝐶𝜕𝑥²−𝑣𝜕𝐶𝜕𝑥−𝜆𝑅𝐶=𝑅𝜕𝐶().𝜕𝑡+𝑓𝑥,𝑡 (7) Subject to the initial and boundary conditions: 𝐶(𝑥,0)=0,𝐶(0,𝑡)=𝑐₀exp(−𝛼𝑡),𝐶_𝑥(∞,𝑡)=0, (8) where 𝐷 is the dispersion coefficient, 𝑣 is the seepage velocity, 𝑅 is the retardation factor, 𝜆 is the radioactive decay constant, 𝑐₀ is the initial concentration, 𝛼 is a positive constant, and 𝑓(𝑥,𝑡) is any source and sink in the system. However, in the case of the groundwater pollution, the function 𝑓(𝑥,𝑡) is always neglected because we assume that there is no source and sink in the system under investigation. Therefore, in our case we will set the function to be zero. In order to include explicitly the medium structure or external field changes with time into the mathematical formulation, we replace the standard derivative operator by the variational order derivative operator as follows: 𝑓()𝜕𝑥,𝑡+𝐷²𝐶𝜕𝑥²−𝑣𝜕𝐶𝜕𝑥−𝜆𝑅𝐶=𝑅𝐷₀^{𝛼(𝑡,𝑥)}(𝐶). (9) Subject to the initial and boundary conditions: 𝐶(𝑥,0)=0,𝐶(0,𝑡)=𝑐₀exp−𝛾𝑡,𝐶_𝑥(∞,𝑡)=0, (10) where 𝛼(𝑡,𝑥) is a continuous function in (0,1] that can be approximated from the field observation. Note that there is no analytical solution available for this new equation; therefore, our next concern is to establish the existence and uniqueness of such solution. This will be done in the next section.

3. Numerical Solution

Numerical methods acquiesce estimated solutions to the overriding equation all the way through the discretization of space and time. Inside the discredited problem domain, the variable internal properties, boundaries, and stresses of the system are approximated. Deterministic, distributed-parameter, and numerical models can relax the rigid idealized conditions of analytical models or lumped-parameter models, and they can therefore be more realistic and flexible for simulating fields’ conditions [19–29]. The finite difference schemes for constant-order time or space fractional diffusion equations have been widely studied [19–29]. Before performing the numerical methods, we assume that (9) has a unique and sufficiently smooth solution. To establish the numerical schemes for the above equation, we let 𝑥_𝑙=𝑙ℎ, 0≤𝑙≤𝑀, 𝑀ℎ=𝐿, 𝑡_𝑘=𝑘𝜏, 0≤𝑘≤𝑁, and 𝑁𝜏=𝑇, and we let ℎ be the step and 𝜏 be the time size, and 𝑀 and 𝑁 be grid points.

3.1. Crank-Nicolson Scheme [30]

We bring in the Crank-Nicolson idea as follows. At the outset, the discretization of first- and second-order space derivative is stated as 𝜕𝐶=1𝜕𝑥2𝐶𝑥_𝑙+1,𝑡_𝑘+1−𝐶𝑥_𝑙−1,𝑡_𝑘+12(ℎ)+𝐶𝑥_𝑙+1,𝑡_𝑘−𝐶𝑥_𝑙−1,𝑡_𝑘2(ℎ)(ℎ),𝜕+𝑂₂𝐶𝜕𝑥²=12𝐶𝑥_𝑙+1,𝑡_𝑘+1−2𝐶𝑥_𝑙,𝑡_𝑘+1+𝐶𝑥_𝑙−1,𝑡_𝑘+1(ℎ)²+𝐶𝑥_𝑙+1,𝑡_𝑘−2𝐶𝑥_𝑙,𝑡_𝑘+𝐶𝑥_𝑙−1,𝑡_𝑘(ℎ)²+𝑂ℎ²,1𝐶=2𝐶𝑥_𝑙,𝑡_𝑘+1+𝐶𝑥_𝑙,𝑡_𝑘. (11)(12)(13) The Crank-Nicolson scheme for the VO time fractional diffusion model can be stated as follows: 𝜕^{𝛼𝑙𝑘+1}𝐶𝑥_𝑙,𝑡_𝑘+1𝜕𝑡^{𝛼𝑙𝑘+1}=𝜏^{−𝛼𝑙𝑘+1}Γ2−𝛼_𝑙^𝑘+1𝐶𝑥_𝑙,𝑡_𝑘+1−𝐶𝑥_𝑙,𝑡_𝑘+^𝑘_𝑗=1𝐶𝑥_𝑙,𝑡_{𝑘+1−𝑗}−𝐶𝑥_𝑙,𝑡_𝑘−𝑗×𝑗+1^{1−𝛼𝑙𝑘+1}−𝑗^{1−𝛼𝑙𝑘+1}+𝑂(𝜏). (14) Now replacing (11), (12), (13), and (14) in (9) we obtain the following: 𝑏_𝑙^𝑘+1𝐶^𝑘+1_𝑙+1−2𝐶_𝑙^𝑘+1+𝐶^𝑘_𝑙−1+𝐶^𝑘_𝑙+1−2𝐶^𝑘_𝑙+𝐶^𝑘_𝑙−1+𝑐_𝑙^𝑘+1𝐶^𝑘+1_𝑙+1−𝐶^𝑘+1_𝑙−1+𝐶^𝑘_𝑙+1−𝐶^𝑘_𝑙−1+𝑑_𝑙^𝑘+1𝐶_𝑙^𝑘+1−𝐶^𝑘_𝑙=𝐶_𝑙^𝑘+1−𝐶^𝑘_𝑙+^𝑘_𝑖=1𝐶_𝑙^{𝑘+1−𝑖}−𝐶_𝑙^𝑘−𝑖𝑒_𝑗^𝑙,𝑘+1+𝜏^{𝛼𝑙𝑘+1}Γ2−𝛼_𝑙^𝑘+1𝑓_𝑙^𝑘+1. (15) Here 𝑏_𝑙^𝑘+1=𝐷𝜏^{𝛼𝑙𝑘+1}Γ2−𝛼_𝑙^𝑘+12𝑅(ℎ)²,𝑐_𝑙^𝑘+1=−𝑣𝜏^{𝛼𝑙𝑘+1}Γ2−𝛼_𝑙^𝑘+1,𝑑2ℎ𝑅_𝑙𝑘+1=𝜆𝜏^{𝛼𝑙𝑘+1}Γ2−𝛼_𝑙^𝑘+12 (16)𝑒_𝑗^𝑙,𝑘+1=(𝑗+1)^{1−𝛼𝑙𝑘+1}−(𝑗)^{1−𝛼𝑙𝑘+1}, 𝑝_𝑗^𝑙,𝑘+1=𝑒^𝑙,𝑘+1_𝑗−1−𝑒_𝑗^𝑙,𝑘+1𝑗=1,2,…,𝑁, and 𝑙=1,2,…,𝑀−1. It is important to point out that the sum term on the right-hand side of (15) automatically vanishes when 𝑘=0. Then, (15) can be reformulated as: 𝑏¹_𝑙𝐶¹_𝑙+1−2𝐶¹_𝑙−1+𝐶⁰_𝑙−1+𝐶⁰_𝑙+1−2𝐶⁰_𝑙−1+𝐶⁰_𝑙−1+𝑐¹_𝑙𝐶¹_𝑙+1−𝐶¹_𝑙−1+𝐶⁰_𝑙+1−𝐶⁰_𝑙−1+𝑑¹_𝑙𝐶¹_𝑙−𝐶⁰_𝑙=𝜏^𝛼1𝑙Γ2−𝛼¹_𝑙𝑓¹_𝑙,𝑏for𝑘=0,_𝑙𝑘+1𝐶^𝑘+1_𝑙+1−2𝐶_𝑙^𝑘+1+𝐶^𝑘_𝑙−1+𝐶^𝑘_𝑙+1−2𝐶^𝑘_𝑙+𝐶^𝑘_𝑙−1+𝑐_𝑙^𝑘+1𝐶^𝑘+1_𝑙+1−𝐶^𝑘+1_𝑙−1+𝐶^𝑘_𝑙+1−𝐶^𝑘_𝑙−1+𝑑_𝑙^𝑘+1𝐶_𝑙^𝑘+1−𝐶^𝑘_𝑙=𝐶_𝑙^𝑘+1−𝐶^𝑘_𝑙+^𝑘_𝑖=1𝐶_𝑙^{𝑘+1−𝑖}−𝐶_𝑙^𝑘−𝑖𝑒_𝑗^𝑙,𝑘+1+𝜏^{𝛼𝑙𝑘+1}Γ2−𝛼_𝑙^𝑘+1𝑓_𝑙^𝑘+1,for𝑘≥1. (17)(18) The above discretization can be rewritten in the following matrix form.

3.2. Stability Analysis of the Crank-Nicolson Scheme

In this section, we will analyze the stability conditions of the Crank-Nicolson scheme for the generalized advection dispersion equation.

Let 𝜁^𝑘_𝑙=𝑐^𝑘_𝑙−𝐶^𝑘_𝑙; here 𝐶^𝑘_𝑙 is the approximate solution at the point (𝑥_𝑙,𝑡_𝑘), (𝑘=1,2…,𝑁,𝑙=1,2,…,𝑀−1) and in addition 𝜁^𝑘=[𝜁^𝑘₁,𝜁^𝑘₂,…,𝜁^𝑘_𝑀−1]^𝑇 and the function 𝜁^𝑘(𝑥) is chosen to be 𝜁^𝑘(𝑥)=𝜁^𝑘_𝑙if𝑥_𝑙−ℎ2<𝑥≤𝑥_𝑙+ℎ2ℎ,𝑙=1,2,…,𝑀−10if𝐿−2<𝑥≤𝐿. (19) Then, the function 𝜁^𝑘(𝑥) can be expressed in Fourier series as follows: 𝜁^𝑘(𝑥)=^𝑚=∞_𝑚=−∞𝛿_𝑚(𝑚)exp2𝑖𝜋𝑚𝑘𝐿,𝛿_𝑘(𝑥)=1𝐿^𝐿₀𝜌^𝑘(𝑥)exp2𝑖𝜋𝑚𝑥𝐿𝑑𝑥. (20) It was established by [25] that 𝜁²²₂=^𝑚=∞_𝑚=−∞𝛿_𝑘(𝑚)². (21) Observe that, for all 𝑘,𝑙≥1, 0≤1−𝛼_𝑙^𝑘+1<1; in addition, according to the problem in point, the velocity seepage 𝑣, the dispersion coefficient 𝐷, the retardation factor 𝑅, and the radioactive decay constant 𝜆 are positive constants. Then the following properties of the coefficients 𝑐_𝑙^𝑘+1, 𝑑_𝑙^𝑘+1, 𝑝_𝑗^𝑙,𝑘+1, and 𝑏_𝑙^𝑘+1 can be established as follows.(1)𝑑_𝑙^𝑘+1, 𝑏_𝑙^𝑘+1 are positive for all 𝑙=1,2…,𝑀−1.(2)𝑐_𝑙^𝑘+1 is negative for all 𝑙=1,2…,𝑀−1.(3)0<𝑒_𝑗^𝑙,𝑘≤𝑒^𝑙,𝑘_𝑗−1≤1 for all 𝑙=1,2…,𝑀−1.(4)0≤𝑝_𝑗^𝑙,𝑘≤1, ∑^𝑘−1_𝑗=0𝑝^𝑙,𝑘+1_𝑗+1=1−𝑒_𝑘^𝑙,𝑘+1 for all 𝑙=1,2…,𝑀−1.

The error committed while approximating the solution of the generalized advection dispersion equation with Crank-Nicolson scheme can be presented as follows: 𝑏_𝑙^𝑘+1𝜁^𝑘+1_𝑙+1−2𝜁_𝑙^𝑘+1+𝜁^𝑘+1_𝑙−1+𝜁^𝑘_𝑙+1−2𝜁^𝑘_𝑙+𝜁^𝑘_𝑙−1+𝑐_𝑙^𝑘+1𝜁^𝑘+1_𝑙+1−𝜁^𝑘+1_𝑙−1+𝜁^𝑘_𝑙+1−𝜁^𝑘_𝑙−1+𝑑_𝑙^𝑘+1𝜁_𝑙^𝑘+1−𝜁^𝑘_𝑙=𝜁_𝑙^𝑘+1−𝜁^𝑘_𝑙+^𝑘_𝑖=1𝜁_𝑙^{𝑘+1−𝑖}−𝜁_𝑙^𝑘−𝑖𝑒_𝑗^𝑙,𝑘+1+𝜏^{𝛼𝑙𝑘+1}Γ2−𝛼_𝑙^𝑘+1𝑓_𝑙^𝑘+1. (22) If we assume that 𝜁^𝑘_𝑙 in (22) can be put in the delta-exponential form as follows: 𝜁^𝑘_𝑙=𝛿_𝑘exp𝑖𝜑𝑙𝑘, (23) where 𝜑 is a real spatial wave number, replacing (23) in (22) we obtain 1+4𝑏¹_𝑙sin²𝜑ℎ2−2𝑑¹_𝑙sin²𝜑ℎ2𝛿₁=1−4𝑏¹_𝑙sin²𝜑ℎ2−2𝑑¹_𝑙sin²𝜑ℎ2𝛿₀,for𝑘=0,1+4𝑏_𝑙^1+𝑘sin²𝜑ℎ2−2𝑑_𝑙^1+𝑘sin²𝜑ℎ2𝛿_𝑘+1=1−4𝑏_𝑙^1+𝑘sin²𝜑ℎ2−2𝑑_𝑙^1+𝑘sin²𝜑ℎ2−𝑒₁^𝑙,𝑘+1𝛿_𝑘+^𝑘−1_𝑗=0𝑝^𝑙,𝑘+1_𝑗+1𝛿_𝑘−𝑗+𝑒_𝑘^𝑙,𝑘+1𝛿₀,for𝑘=1,2…,𝑁−1. (24) Equation (24) can be written in the following form: 𝛿₁=1−4𝑏¹_𝑙sin²𝜑ℎ/2−2𝑑¹_𝑙sin²𝜑ℎ/2𝛿₀1+4𝑏¹_𝑙sin²𝜑ℎ/2−2𝑑¹_𝑙sin²,𝛿𝜑ℎ/2_𝑘+1=1−4𝑏_𝑙^1+𝑘sin²𝜑ℎ2−2𝑑_𝑙^1+𝑘sin²𝜑ℎ2−𝑒₁^𝑙,𝑘+1𝛿_𝑘+^𝑘−1_𝑗=0𝑝^𝑙,𝑘+1_𝑗+1𝛿_𝑘−𝑗+𝑒_𝑘^𝑙,𝑘+1𝛿₀×1+4𝑏_𝑙^1+𝑘sin²𝜑ℎ2−2𝑑_𝑙^1+𝑘sin²𝜑ℎ2⁻¹. (25) Our next concern here is to show that for all 𝑘=1,2,…,𝑁−1 the solution of (25) satisfies the following condition: 𝛿_𝑘<𝛿₀. (26) To achieve this, we make use of the recurrence technique on the natural number 𝑘.

For 𝑘=1 and remembering that 𝑑_𝑙^𝑘+1, 𝑏_𝑙^𝑘+1 are positive for all 𝑙=1,2…,𝑀−1, then we obtain 𝛿₁𝛿₀=1−4𝑏¹_𝑙sin²𝜑ℎ/2−2𝑑¹_𝑙sin²𝜑ℎ/21+4𝑏¹_𝑙sin²𝜑ℎ/2−2𝑑¹_𝑙sin²𝜑ℎ/2<1. (27) Assuming that for 𝑚=2,3,…,𝑘 the property is verified, then 𝛿_𝑘+1=1−4𝑏_𝑙^1+𝑘sin²𝜑ℎ2−2𝑑_𝑙^1+𝑘sin²𝜑ℎ2−𝑒₁^𝑙,𝑘+1𝛿_𝑘+^𝑘−1_𝑗=0𝑝^𝑙,𝑘+1_𝑗+1𝛿_𝑘−𝑗+𝑒_𝑘^𝑙,𝑘+1𝛿₀×1+4𝑏_𝑙^1+𝑘sin²𝜑ℎ2−2𝑑_𝑙^1+𝑘sin²𝜑ℎ2⁻¹. (28) Making use of the triangular inequality, we obtain 𝛿_𝑘+1≤1−4𝑏_𝑙^1+𝑘sin²𝜑ℎ2−2𝑑_𝑙^1+𝑘sin²𝜑ℎ2−𝑒₁^𝑙,𝑘+1𝛿_𝑘+^𝑘−1_𝑗=0𝑝^𝑙,𝑘+1_𝑗+1𝛿_𝑘−𝑗+𝑒_𝑘^𝑙,𝑘+1𝛿₀×1+4𝑏_𝑙^1+𝑘sin²𝜑ℎ2−2𝑑_𝑙^1+𝑘sin²𝜑ℎ2⁻¹. (29) Using the recurrence hypothesis, we have 𝛿_𝑘+1≤1−4𝑏_𝑙^1+𝑘sin²𝜑ℎ2−2𝑑_𝑙^1+𝑘sin²𝜑ℎ2+^𝑘−1_𝑗=0𝑝^𝑙,𝑘+1_𝑗+1×1+4𝑏_𝑙^1+𝑘sin²𝜑ℎ2−2𝑑_𝑙^1+𝑘sin²𝜑ℎ2⁻¹𝛿₀,𝛿_𝑘+1≤1+4𝑏_𝑙^1+𝑘sin²𝜑ℎ/2−2𝑑_𝑙^1+𝑘sin²𝜑ℎ/21+4𝑏_𝑙^1+𝑘sin²𝜑ℎ/2−2𝑑_𝑙^1+𝑘sin²𝛿𝜑ℎ/2⁰,𝛿_𝑘+1≤𝛿₀ (30) and this completes the proof.

3.3. Convergence Analysis of the Crank-Nicolson Scheme

If we assume that 𝑐(𝑥_𝑙,𝑡_𝑘)(𝑙=1,2,…,𝑀,𝑘=1,2…,𝑁−1) is the exact solution of our problem at the point (𝑥_𝑙,𝑡_𝑘), by letting 𝑤^𝑘_𝑙=𝑐(𝑥_𝑙,𝑡_𝑘)−𝐶^𝑘_𝑙 and 𝑤^𝑘=(0,𝑤^𝑘₁,𝑤^𝑘₂,…,𝑤^𝑘_𝑀−1) and substituting this in (18), we obtain 𝑏¹_𝑙𝑤¹_𝑙+1−2𝑤¹_𝑙−1+𝑐¹_𝑙𝑤¹_𝑙+1−𝑤𝐶¹_𝑙−1+𝑑¹_𝑙𝑤¹_𝑙=𝐹¹_𝑙,𝑏for𝑘=0,_𝑙𝑘+1𝑤^𝑘+1_𝑙+1−2𝑤_𝑙^𝑘+1+𝑤^𝑘+1_𝑙−1+𝑐_𝑙^𝑘+1𝑤^𝑘+1_𝑙+1−𝑤^𝑘+1_𝑙−1+𝑑_𝑙^𝑘+1𝑤_𝑙^𝑘+1=𝐹_𝑙^𝑘+1+^𝑘_𝑖=1𝑤_𝑙^𝑘−𝑖𝑒_𝑗^𝑙,𝑘+1,for𝑘≥1. (31) Here 𝐹_𝑙^𝑘+1=𝑑_𝑙^𝑘+1−𝐶𝑥_𝑙,𝑡_𝑘+1+𝐶𝑥_𝑙,𝑡_𝑘+^𝑘_𝑖=1𝐶𝑥_𝑙,𝑡_{𝑘−𝑖+1}𝑒_𝑗^𝑙,𝑘+1+𝑏_𝑙^𝑘+1−𝐶𝑥_𝑙+1,𝑡_𝑘+2𝐶𝑥_𝑙,𝑡_𝑘−𝐶𝑥_𝑙−1,𝑡_𝑘−𝑐_𝑙^𝑘+1𝐶𝑥_𝑙+1,𝑡_𝑘−𝐶𝑥_𝑙−1,𝑡_𝑘. (32) From (11), (12), and (14), 𝜕𝐶𝑥_𝑙,𝑡_𝑘+1𝜕𝑥+𝑉₁=12𝐶𝑥_𝑙+1,𝑡_𝑘+1−𝐶𝑥_𝑙−1,𝑡_𝑘+12(ℎ)+𝐶𝑥_𝑙+1,𝑡_𝑘−𝐶𝑥_𝑙−1,𝑡_𝑘2(ℎ),𝜕²𝐶𝑥_𝑙,𝑡_𝑘+1𝜕𝑥²+ℎ²𝑉₂=12𝐶𝑥_𝑙+1,𝑡_𝑘+1−2𝐶𝑥_𝑙,𝑡_𝑘+1+𝐶𝑥_𝑙−1,𝑡_𝑘+1ℎ²+𝐶𝑥_𝑙+1,𝑡_𝑘−2𝐶𝑥_𝑙,𝑡_𝑘+𝐶𝑥_𝑙−1,𝑡_𝑘ℎ²,𝜕^{𝛼𝑙𝑘+1}𝐶𝑥_𝑙,𝑡_𝑘+1𝜕𝑡^{𝛼𝑙𝑘+1}+𝜏𝑉₃=𝜏^{−𝛼𝑙𝑘+1}Γ2−𝛼_𝑙^𝑘+1𝐶𝑥_𝑙,𝑡_𝑘+1−𝐶𝑥_𝑙,𝑡_𝑘+^𝑘_𝑗=1𝐶𝑥_𝑙,𝑡_{𝑘+1−𝑗}−𝐶𝑥_𝑙,𝑡_𝑘−𝑗×𝑗+1^{1−𝛼𝑙𝑘+1}−𝑗^{1−𝛼𝑙𝑘+1}. (33) From the above we have that 𝐹_𝑙^𝑘+1≤𝑉𝜏^{1+𝛼𝑙𝑘+1}+ℎ²𝜏^𝛼𝑘𝑙+ℎ𝜏^𝛼𝑘𝑙, (34) where 𝑉₁, 𝑉₂, 𝑉₃, and 𝑉 are constants. Taking into account Caputo type fractional derivative, the detailed error analysis on the above schemes can refer to the work by Diethelm et al. [31] and further work by Li and Tao [32].

Lemma 4. ‖𝑤^𝑘+1‖_∞≤𝑉(𝜏^{1+𝛼𝑙𝑘+1}+ℎ²𝜏^𝛼𝑘𝑙+ℎ𝜏^𝛼𝑘𝑙)(𝑒_𝑗^𝑙,𝑘+1)⁻¹ is true for (𝑘=0,1,2…,𝑁−1) where ‖𝑤^𝑘‖_∞=max_{1≤𝑙≤𝑀−1}(𝑤^𝑘) and 𝑉 is a constant. In addition, 𝛼^𝑘+1=min_{1≤𝑙≤𝑀−1}𝛼_𝑙^𝑘+1,𝑖𝑓𝜏<1max_{1≤𝑙≤𝑀−1}𝛼_𝑙^𝑘+1,𝑖𝑓𝜏>1. (35) This can be achieved via the recurrence technique on the natural number 𝑘.
When 𝑘=0, one has the following: 𝑤¹_𝑙≤𝑐¹_𝑙+𝑏¹_𝑙𝑤¹_𝑙+1+𝑐¹_𝑙−2𝑏¹_𝑙𝑤¹_𝑙−1+𝑑¹_𝑙𝑤¹_𝑙=𝐹¹_𝑙≤𝑉𝜏^{1+𝛼𝑙𝑘+1}+ℎ²𝜏^𝛼𝑘𝑙+ℎ𝜏^𝛼𝑘𝑙𝑒_𝑗^{𝑙,𝑘+1−1}. (36) Now suppose that ‖𝑤^𝑖+1‖_∞≤𝑉(𝜏^{1+𝛼𝑙𝑖+1}+ℎ²𝜏^𝛼𝑖𝑙+ℎ𝜏^𝛼𝑖𝑙)(𝑒_𝑗^𝑙,𝑖+1)⁻¹, 𝑖=1,…,𝑁−2. Then 𝑤_𝑙^1+𝑘≤𝑏_𝑙^𝑘+1𝑤^𝑘+1_𝑙+1−2𝑤_𝑙^𝑘+1+𝑤^𝑘+1_𝑙−1+𝑐_𝑙^𝑘+1𝑤^𝑘+1_𝑙+1−𝑤^𝑘+1_𝑙−1+𝑑_𝑙^𝑘+1𝑤_𝑙^𝑘+1≤𝑏_𝑙^𝑘+1+𝑐_𝑙^𝑘+1𝑤^𝑘+1_𝑙+1+𝑏_𝑙^𝑘+1−𝑐_𝑙^𝑘+1𝑤^𝑘+1_𝑙−1+𝑑_𝑙^𝑘+1−2𝑏_𝑙^𝑘+1𝑤_𝑙^𝑘+1=𝐹_𝑙^𝑘+1+^𝑘_𝑖=1𝑤_𝑙^𝑘−𝑖𝑒_𝑗^𝑙,𝑘+1≤𝐹_𝑙^𝑘+1+^𝑘_𝑖=1𝑤_𝑙^𝑘−𝑖𝑒_𝑗^𝑙,𝑘+1≤𝑉𝜏^{1+𝛼𝑙𝑘+1}+ℎ²𝜏^𝛼𝑘𝑙+ℎ𝜏^𝛼𝑘𝑙+^𝑘_𝑖=1𝑤_𝑙^𝑘−𝑖_∞𝑒_𝑗^𝑙,𝑘+1≤𝑉𝜏^{1+𝛼𝑙𝑘+1}+ℎ²𝜏^𝛼𝑘𝑙+ℎ𝜏^𝛼𝑘𝑙×𝑒_𝑗^𝑙,𝑘+1+𝑒₀^𝑙,𝑘+1−𝑒_𝑗^𝑙,𝑘+1𝑒_𝑗^{𝑙,𝑘+1−1}≤𝑉𝜏^{1+𝛼𝑙𝑘+1}+ℎ²𝜏^𝛼𝑘𝑙+ℎ𝜏^𝛼𝑘𝑙𝑒₀^𝑙,𝑘+1𝑒_𝑗^{𝑙,𝑘+1−1},𝑤_𝑙^1+𝑘≤𝑉𝜏^{1+𝛼𝑙𝑘+1}+ℎ²𝜏^𝛼𝑘𝑙+ℎ𝜏^𝛼𝑘𝑙𝑒_𝑗^{𝑙,𝑘+1−1}. (37)This completes the proof.

Theorem 5. The Crank-Nicolson scheme is convergent, and there exists a positive constant 𝑉 such that 𝐶^𝑘_𝑙−𝑐𝑥_𝑙,𝑡_𝑘≤𝑉𝜏+ℎ+ℎ²,𝑙=1,2,…,𝑀−1,𝑘=1,2,…,𝑁. (38)

An interested reader can find the solvability of the Crank-Nicholson scheme in the work done by [33]. Therefore, the details of the proof will not be presented in this paper.

4. Numerical Simulations

In this section, we present the numerical solutions of the following three cases.

The case where we assume that the variation of the geological formation in which the contamination is migrating can be modelled as function of time and space as 𝛼(𝑥,𝑡). To include explicitly this variability in the mathematical formulation, we reformulate the hydrodynamic dispersion equation as follows: 𝐷𝜕²𝐶𝜕𝑥²−𝑣𝜕𝐶𝜕𝑥−𝜆𝑅𝐶=𝑅𝐷₀^{𝛼(𝑡,𝑥)}(𝐶),𝛼()(),𝐶()()𝑥,𝑡=1−sin𝑥𝑡𝑥,0=0,𝐶0,𝑡=𝑐₀exp−𝛾𝑡. (39) The theoretical parameters used in this simulation are given in Table 1.

Table 1 

Theoretical parameters.

Figure 1 shows the approximate solution of the main problem as function of space and time. Figure 2 shows the density plot of the approximate solution of the main problem. Figure 3 shows the variation of the concentration as function of time for a fixed space. Figures 1 and 3 show the solution of the variation order advection and advection equation, respectively, and Figure 2 shows the contour map of variational order advection equation.

Figure 1 

Approximate solution of the modified equation as function of time and space.

Figure 2 

Solution of the modified equation seing as contour map.

Figure 3 

Solution of the standard equation.

5. Discussion and Conclusions

Natural geological deposits with highly contrasting permeability may form mobile and relatively immobile zones, where the potential mass exchange between mobile and immobile zones results in a wide time distribution for solute “trapping.” The transport process, groundwater, is, by its very nature, always in contact with the matrix of an aquifer. There is thus a possibility that the solutes may interact with the rock matrix and one another. A true mathematical model for groundwater pollution must therefore be able to account for interactions between the dissolved solids and matrix of the aquifer. It will thus be advantageous to look at the nature of the interactions between dissolved solids and a porous medium that may be expected in groundwater pollution. Experimental evidence indicates that when a dissolved solid comes in contact with the matrix of a porous medium it may (a) pass through the medium with no apparent effect, (b) be absorbed by the porous matrix, and (c) react with the porous matrix and other substances dissolved in the fluid. The dissolved solids encountered in porous flow are, for this reason, often classified as conservative, nonconservative, and reactive tracers. This behaviour implies that the quantity of dissolved solids in a porous medium depends not only on the flow pattern but also on the nature of the porous matrix and the solution. These situations (a), (b), and (c) can be characterized efficiently by the time-nonlocal model, including the time variational order advection dispersion equation. If the high-permeable material tends to form preferential flow paths, such as the interconnected paleochannels observed in alluvial depositional systems, then the solute transport may show a heavy leading edge, which can be described by the Variational Order Advection Dispersion Equation (VOADE). Development of partial differential equations such as the advection-dispersion equation (ADE) begins with assumptions about the random behaviour of a single particle: possible velocities it may experience in a flow field and the length of time it may be immobilized. When assumptions underlying the ADE are relaxed, a fractional ADE (FADE) can arise, with a non-integer-order derivative on time or space terms. Fractional ADEs are nonlocal; they describe transport affected by hydraulic conditions at a distance. Space fractional ADEs arise when velocity variations are heavy tailed and describe particle motion that accounts for variation in the flow field over the entire system [34]. But when the geological formation, under which this pollution moves, is changing in time and space, the VOADE can be used. Many other uses of the variable-order derivatives can be found in the following works [35–40].

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgement

This work was partially financially supported by the Leon Claude Fellowships of 2014 of South Africa. Adem Kilicman gratefully acknowledges that this research was partially supported by the University Putra Malaysia under the Grant Scheme GB-IBT having the project number GB-IBT/2013/9420100.