Stability and Convergence of a Time-Fractional Variable Order Hantush Equation for a Deformable Aquifer

and Applied Analysis 3 3. Numerical Solution Environmental phenomena, such as groundwater flow described by variational order derivative, are highly complex phenomena, which do not lend themselves readily to analysis of analytical models. The discussion presented in this section will therefore be devoted to the derivation of numerical solution to the modified Hantush equation (6). Solving difficult equations with numerical scheme has been passionate exercise for many scholars [9–20]. However, there exists numerous of this scheme in the literature [14– 20]. Some of these numerical techniques are very accurate while approximating solutions of difficult equations. These numerical methods yield approximate solutions to the governing equation through the discretization of space and time [7]. Within the discredited problem domain, the variable internal properties, boundaries, and stresses of the system are approximated. Deterministic, distributed-parameter, numerical models can relax the rigid idealised conditions of analytical models or lumped-parameter models, and they can therefore be more realistic and flexible for simulating fields conditions [7]. Recently Atangana and Botha [7] have extended the groundwater flow model to the concept of time-fractional variable order derivative; they solved the resulting equation via Crank-Nicolson numerical scheme. The finite difference schemes for constant-order timeor space-fractional diffusion equations have beenwidely studied in [9–14]. Recently, Sun et al. [21] studied the solution of the advection dispersion equation with time-fractional variable order derivative. The study of the implicit difference approximation scheme for constant-order time-fractional diffusion equations was presented in [15]. Recently, the weighted average finite difference method was introduced [16]. The matrix approach for fractional diffusion equations was proposed [17], and Hanert proposed a flexible numerical scheme for the discretization of the space-time fractional diffusion equation (see [18]). Recently, the numerical scheme for VO space-fractional advection-dispersion equation was considered [19]. The investigation of the explicit scheme for VO nonlinear space-fractional diffusion equation was done (see [20]). 3.1. Crank-Nicolson Scheme [22]. Before performing the numerical methods, we assume that (3) has a unique and sufficiently smooth solution. To establish the numerical schemes for the above equation, we let x l = lh, 0 ≤ l ≤ M, Mh = L, t k = kτ, 0 ≤ k ≤ N, Nτ = T, h is the step and τ is the time size, andM andN are grid points. We introduce the Crank-Nicolson scheme as follows. Firstly, the discretization of firstand second-order space derivative is stated as ∂s ∂r = 1 2 (( s (r l+1 , t k+1 ) − Φ (r l−1 , t k+1 ) 2 (h) ) +( s (r l+1 , t k ) − s (r l−1 , t k ) 2 (h) )) + O (h) , (7) ∂2s ∂r = 1 2 (( s (r l+1 , t k+1 ) − 2s (r l , t k+1 ) + s (r l−1 , t k+1 ) (h)2 ) +( s (r l+1 , t k ) − 2s (r l , t k ) + s (r l−1 , t k ) (h)2 ))


Introduction
An aquifer is an underground layer of water-bearing permeable rock or unconsolidated materials (gravel, sand, or silt) from which groundwater can be extracted using a water well.The study of water flow in aquifers and the characterization of aquifers is called hydrogeology.Related terms include aquitard, which is a bed of low permeability along an aquifer, see [1] and aquiclude (or aquifuge), which is a solid, impermeable area underlying or overlying an aquifer.If the impermeable area overlies the aquifer, pressure could cause it to become a confined aquifer.There are two end members in the spectrum of types of aquifers: confined and unconfined (with semiconfined being in between).Unconfined aquifers are sometimes also called water table or phreatic aquifers, because their upper boundary is the water table or phreatic surface.Typically but not always, the shallowest aquifer at a given location is unconfined, meaning that it does not have a confining layer (an aquitard or an aquiclude) between it and the surface.When a leaky aquifer is pumped, the piezometric level of the aquifer in the well is lowered.This lowering spreads radially outward as pumping continues, creating a difference in hydraulic head between the aquifer and the aquitards.Consequently, the groundwater in the aquitards will start moving vertically downward to join the water in the aquifer.The aquifer is thus partially recharged by downward percolation from the aquitards.As pumping continues, the percentage of the total discharge derived from this percolation increases.After a certain period of pumping, equilibrium will be established between the discharge rate of the pump and the recharge rate by vertical flow through the aquitards.This steady state will be maintained as long as the water table in the aquitards is kept constant.Figure 1 shows the piezometric level after the start of pumping in a leaky aquifer.
Hantush was the first to derive a partial differential equation describing such phenomena.However, due to the deformation of some aquifers, the Hantush equation is not able to account for the effect of the changes in the mathematical formulation.The purpose of this work is therefore devoted to the discussion underpinning the description of the groundwater flow through the deformable leaky aquifer, on one hand.On the other hand, we present the derivation of the approximate solution of the modified equation via the Crank-Nicolson scheme.We will start with the definition of the variational order derivative and the problem modification.

Definition and Problem Modification
For the readers that are not acquainted with the concept of the variational order derivative, we start this section.We present the basic definition of this derivative.

Variational Order of Differential
Operator.Let  : R → R,  → () denote a continuous but necessary differentiable; let () be a continuous function in (0, 1].Then its variational order differential is defined as The above derivative is called the Caputo variational order differential operator; additionally the derivative of the constant is zero.

Problem Formulation.
Groundwater models describe the groundwater flow and transport processes using mathematical equations based on certain simplifying assumptions.These assumptions typically involve the direction of flow, geometry of the aquifer, and the heterogeneity or anisotropy of sediments or bedrock within the aquifer.This geological formation through which the groundwater flows changes in time and space.
The simplest generalization of groundwater flow equation, which incidentally is also in accord with true physics of the phenomenon, is to assume that water level is not in a steady but transient state.In 1935, Theis [2] was the first to develop a formula for unsteady-state flow that introduces the time factor and the storativity.He noted that when a well penetrating an extensive confined aquifer is pumped at a constant rate, the influence of the discharge extends outward with time.The rate of decline of head, multiplied by the storativity and summed over the area of influence, equals the discharge.The unsteady-state (or Theis) equation, which was derived from the analogy between the flow of groundwater and the conduction of heat, is perhaps the most widely used partial differential equation in groundwater investigations The above equation is classified under parabolic equations.However, very few geological formations are completely impermeable to fluids.Leakage of the water could thus occur, should a confined aquifer be over-or underlain by another aquifer.The behaviour of such an aquifer, often referred to as a leaky or semiconfined aquifer, needs thus not be the same as that of a confined aquifer.Although the nature of a semiaquifer differs from that of a true aquifer, it is still possible to use the basic principles of confined flow to arrive at the governing equation for such aquifer.This is in particular true in those situations where the confining layer between the two aquifers is not too thick and the flow is mainly in the vertical direction.
According to Hantush and Jacob [3,4], the drawdown due to pumping a leaky aquifer can be described by the following equation: where  is the drawdown or change in the level of water;  is the specific storativity of the aquifer, and  is the transmissivity: with  and   as the hydraulic conductivities of the main and confining layers, respectively,  and  are the thicknesses of the main and confining layers, respectively, and  is the discharge rate of the pumping.This partial differential equation describing the movement of water through the geological formation during the pumping is subjected to the following initial and boundary conditions: However, when we consider the diffusion process in the porous medium, if the medium structure or external field changes with time, in this situation, the ordinary integerorder and constant-order fractional diffusion equation model cannot be used to well characterize such phenomenon [5,6].This is the case of the groundwater flow in the deformable aquifer, the medium through which the flow occurs changes with time and space [7,8].Feature that equation Hantush cannot handle this case.One of the purposes of this work is therefore devoted to the discussion underpinning the description of water flowing through a deformable leaky aquifer, on one hand.In order to include explicitly the variability of the medium through which the flow takes place, the standard version of the partial derivative with respect to time is replaced here with variable-order (VO) fractional to obtain

Numerical Solution
Environmental phenomena, such as groundwater flow described by variational order derivative, are highly complex phenomena, which do not lend themselves readily to analysis of analytical models.The discussion presented in this section will therefore be devoted to the derivation of numerical solution to the modified Hantush equation ( 6).Solving difficult equations with numerical scheme has been passionate exercise for many scholars [9][10][11][12][13][14][15][16][17][18][19][20].However, there exists numerous of this scheme in the literature [14][15][16][17][18][19][20].Some of these numerical techniques are very accurate while approximating solutions of difficult equations.These numerical methods yield approximate solutions to the governing equation through the discretization of space and time [7].Within the discredited problem domain, the variable internal properties, boundaries, and stresses of the system are approximated.Deterministic, distributed-parameter, numerical models can relax the rigid idealised conditions of analytical models or lumped-parameter models, and they can therefore be more realistic and flexible for simulating fields conditions [7].Recently Atangana and Botha [7] have extended the groundwater flow model to the concept of time-fractional variable order derivative; they solved the resulting equation via Crank-Nicolson numerical scheme.The finite difference schemes for constant-order time-or space-fractional diffusion equations have been widely studied in [9][10][11][12][13][14].Recently, Sun et al. [21] studied the solution of the advection dispersion equation with time-fractional variable order derivative.The study of the implicit difference approximation scheme for constant-order time-fractional diffusion equations was presented in [15].Recently, the weighted average finite difference method was introduced [16].The matrix approach for fractional diffusion equations was proposed [17], and Hanert proposed a flexible numerical scheme for the discretization of the space-time fractional diffusion equation (see [18]).Recently, the numerical scheme for VO space-fractional advection-dispersion equation was considered [19].The investigation of the explicit scheme for VO nonlinear space-fractional diffusion equation was done (see [20]).[22].Before performing the numerical methods, we assume that (3) has a unique and sufficiently smooth solution.To establish the numerical schemes for the above equation, we let   = ℎ, 0 ≤  ≤ , ℎ = ,   = , 0 ≤  ≤ ,  = , ℎ is the step and  is the time size, and  and  are grid points.We introduce the Crank-Nicolson scheme as follows.Firstly, the discretization of first-and second-order space derivative is stated as
It is customary in groundwater investigations to choose a point on the centreline of the pumped borehole as a reference for the observations, and therefore, neither the drawdown nor its derivatives will vanish at the origin, as required [7].In such situation, the distribution of the piezometric head in the aquifer is a decreasing function of the distance from the borehole, the expression (1/  ) → 0 [7].Under this situation, the error committed while approximating the solution of the generalized advection dispersion equation with Crank-Nicolson scheme can be presented as follows: If we assume that Υ   in ( 17) can be put in the deltaexponential form as follows [21]: where  is a real spatial wave number, now replacing the above equation ( 18) in ( 17) we obtain Equation ( 19) can be written in the following form: Our next concern here is to show that for all  = 1, 2, . . .,  − 1 the solution of ( 19) satisfies the following condition: To achieve this, we make use of the recurrence technique on the natural number .
For  = 1, and remembering that  Assuming that for  = 2, 3, . . .,  the property is verified, then Making use of the triangular inequality, we obtain which completes the proof that starts at (21) and ends at (25).
From ( 8) and ( 9), we have From the above, we have that where  1 ,  2 , and  are constants.Taking into account Caputo-type fractional derivative, the detailed error analysis on the above schemes can refer to the work done in [21] and further work done in [21,23].
When  = 0, we have the following: which completes the proof of Lemma 1.
An interesting and detailed research can be found on the solvability of the Crank-Nicolson scheme in the work [7,8,22].Therefore, the details of the proof of Theorem 2 (33) will not be presented in this paper.

Numerical Simulations
In this section, we present the numerical simulation of the solution of the modified Hantush equation obtained via the Crank-Nicolson scheme.Here let us consider the following equation:    hydraulic conductivity of the main aquifer of  = 8 feet/day, and a hydraulic conductivity of the leaky   = 0.4 feet/day.A storativity  = 0.001, the leaky factor is 0.00081/day, and finally the flow rate = 4200 feet/day.The red shows the simulation for a distance of  = 400 feet.The black shows the simulation for a distance of  = 200 feet, and finally the blue is the simulation for  = 3650 feet.Figure 3 shows the numerical simulation of the water flowing in the variable leaky aquifer.

Conclusion
The main purpose of this paper was to consider the deformation of the leaky aquifer in the mathematical formulation.However, there are some leaky aquifers that change in time and space.Feature that the Hantush equation cannot be used to describe such situation.Recently, the variational order derivative was found very useful to describe efficiently such situation.We then modified the Hantush equation by replacing the partial derivative with respect to time by the variational order derivative.The result equation was solved numerically using the Crank-Nicolson scheme.The stability and convergence were presented in details.We compared the numerical simulation together with the observed drawdown from field observation.The comparison shows that the modified equation predicts more accurately the real field observation.

Figure 2
Figure 2 shows the numerical simulation for an aquifer thickness of 1000 feet, a transmissivity  = 10000 feet/day, a

Figure 2 :
Figure 2: Numerical comparison of the observed data and solution of the modified equation for different values of time.

Figure 3 :
Figure 3: Numerical simulation of the water flowing in the leaky aquifer.