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

Authors: Abdon Atangana; S. C. Oukouomi Noutchie
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

The medium through which the groundwater moves varies in time and space. The Hantush equation describes the movement of groundwater through a leaky aquifer. To include explicitly the deformation of the leaky aquifer into the mathematical formulation, we modify the equation by replacing the partial derivative with respect to time by the time-fractional variable order derivative. The modified equation is solved numerically via the Crank-Nicolson scheme. The stability and the convergence in this case are presented in details.


---

## Body

## 1. 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.Figure 1
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.


## 2. 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.


### 2.1. Variational Order of Differential Operator

Letf:ℝ→ℝ,x→f(x) denote a continuous but necessary differentiable; let α(x) be a continuous function in (0, 1]. Then its variational order differential is defined as
(1)D0α(x)(f(x))=1Γ(1-α(x))∫0x(x-t)-α(t)df(t)dtdt.The above derivative is called the Caputo variational order differential operator; additionally the derivative of the constant is zero.


### 2.2. 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
(2)SDtΦ(r,t)=TDrrΦ(r,t)+1rDrΦ(r,t).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:
(3)SDts(r,t)=TDrrs(r,t)+1rDrs(r,t)+s(r,t)λ2,
where s is the drawdown or change in the level of water; S is the specific storativity of the aquifer, and T is the transmissivity:
(4)λ2=BK′,
with K and K′ as the hydraulic conductivities of the main and confining layers, respectively, B and d are the thicknesses of the main and confining layers, respectively, and Q 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:
(5)s(r,0)=Φ0,limr→∞Φ(r,t)=Φ0,Q=2πn/2Γ(n/2)rbn-1Kd3-n∂rs(rb,t).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 integer-order 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
(6)SDtα(x,t)s(r,t)=TDrrs(r,t)+1rDrs(r,t)+s(r,t)λ2,0<α(x,t)≤1.


## 2.1. Variational Order of Differential Operator

Letf:ℝ→ℝ,x→f(x) denote a continuous but necessary differentiable; let α(x) be a continuous function in (0, 1]. Then its variational order differential is defined as
(1)D0α(x)(f(x))=1Γ(1-α(x))∫0x(x-t)-α(t)df(t)dtdt.The above derivative is called the Caputo variational order differential operator; additionally the derivative of the constant is zero.


## 2.2. 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
(2)SDtΦ(r,t)=TDrrΦ(r,t)+1rDrΦ(r,t).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:
(3)SDts(r,t)=TDrrs(r,t)+1rDrs(r,t)+s(r,t)λ2,
where s is the drawdown or change in the level of water; S is the specific storativity of the aquifer, and T is the transmissivity:
(4)λ2=BK′,
with K and K′ as the hydraulic conductivities of the main and confining layers, respectively, B and d are the thicknesses of the main and confining layers, respectively, and Q 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:
(5)s(r,0)=Φ0,limr→∞Φ(r,t)=Φ0,Q=2πn/2Γ(n/2)rbn-1Kd3-n∂rs(rb,t).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 integer-order 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
(6)SDtα(x,t)s(r,t)=TDrrs(r,t)+1rDrs(r,t)+s(r,t)λ2,0<α(x,t)≤1.


## 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 time- or space-fractional diffusion equations have been widely 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 xl=lh,  0≤l≤M,Mh=L,tk=kτ,0≤k≤N,Nτ=T,h  is the step and  τ  is the time size, and M and N are grid points. We introduce the Crank-Nicolson scheme as follows. Firstly, the discretization of first- and second-order space derivative is stated as
(7)∂s∂r=12((s(rl+1,tk+1)-Φ(rl-1,tk+1)2(h))+(s(rl+1,tk)-s(rl-1,tk)2(h)))+O(h),(8)∂2s∂r2=12((s(rl+1,tk+1)-2s(rl,tk+1)+s(rl-1,tk+1)(h)2)+(s(rl+1,tk)-2s(rl,tk)+s(rl-1,tk)(h)2))+O(h2),(9)s=12(s(rl,tk+1)+s(rl,tk)).The Crank-Nicolson scheme for the VO time-fractional diffusion model can be stated as follows:(10)∂αlk+1s(rl,tk+1)∂tαlk+1=τ-αlk+1Γ(2-αlk+1)([(j+1)1-αlk+1-(j)1-αlk+1]s(rl,tk+1)-s(rl,tk)+∑j=1k[s(rl,tk+1-j)-s(rl,tk-j)]×[(j+1)1-αlk+1-(j)1-αlk+1]).Now replacing (7), (8), (9), and (10) in (6), we obtain the following:
(11)[Sτ-αlk+1Γ(2-αlk+1)([(j+1)1-αlk+1-(j)1-αlk+1]s(rl,tk+1)-s(rl,tk)+∑j=1k[s(rl,tk+1-j)-s(rl,tk-j)]×[(j+1)1-αlk+1-(j)1-αlk+1])Sτ-αlk+1Γ(2-αlk+1)]=T[12((s(rl+1,tk+1)-2s(rl,tk+1)+s(rl-1,tk+1)(h)2)+(s(rl+1,tk)-2s(rl,tk)+s(rl-1,tk)(h)2))]+1rl[+(s(rl+1,tk)-s(rl-1,tk)2(h)))12((s(rl+1,tk+1)-s(rl-1,tk+1)2(h))+(s(rl+1,tk)-s(rl-1,tk)2(h)))]+12λ2(s(rl,tk+1)+s(rl,tk)).For simplicity, let us put(12)bjl,k+1=(j+1)1-αlk+1-(j)1-αlk+1,Tlk+1=Γ(2-αlk+1)ταlk+1Sh2T,Glk+1=Γ(2-αlk+1)ταlk+1Sh,λjl,k+1=bj-1l,k+1-bjl,k+1,Γlk+1=Γ(2-αlk+1)ταlk+12λ2.Equation (11) becomes
(13)slk+1(1+2Tlk+1-Γlk+1)=sl+1k+1(Tlk+1+Glk+1rl)+sl-1k+1(Tlk+1-Glk+1rl)+sl+1k(Tlk+1-Glk+1rl)+slk+1(1+2Tlk+1)+∑j=1k(slk+1-j-slk-j)λjl,k+1Glk+1+slk(1-Γlk+1+Tlk+1).


## 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 xl=lh,  0≤l≤M,Mh=L,tk=kτ,0≤k≤N,Nτ=T,h  is the step and  τ  is the time size, and M and N are grid points. We introduce the Crank-Nicolson scheme as follows. Firstly, the discretization of first- and second-order space derivative is stated as
(7)∂s∂r=12((s(rl+1,tk+1)-Φ(rl-1,tk+1)2(h))+(s(rl+1,tk)-s(rl-1,tk)2(h)))+O(h),(8)∂2s∂r2=12((s(rl+1,tk+1)-2s(rl,tk+1)+s(rl-1,tk+1)(h)2)+(s(rl+1,tk)-2s(rl,tk)+s(rl-1,tk)(h)2))+O(h2),(9)s=12(s(rl,tk+1)+s(rl,tk)).The Crank-Nicolson scheme for the VO time-fractional diffusion model can be stated as follows:(10)∂αlk+1s(rl,tk+1)∂tαlk+1=τ-αlk+1Γ(2-αlk+1)([(j+1)1-αlk+1-(j)1-αlk+1]s(rl,tk+1)-s(rl,tk)+∑j=1k[s(rl,tk+1-j)-s(rl,tk-j)]×[(j+1)1-αlk+1-(j)1-αlk+1]).Now replacing (7), (8), (9), and (10) in (6), we obtain the following:
(11)[Sτ-αlk+1Γ(2-αlk+1)([(j+1)1-αlk+1-(j)1-αlk+1]s(rl,tk+1)-s(rl,tk)+∑j=1k[s(rl,tk+1-j)-s(rl,tk-j)]×[(j+1)1-αlk+1-(j)1-αlk+1])Sτ-αlk+1Γ(2-αlk+1)]=T[12((s(rl+1,tk+1)-2s(rl,tk+1)+s(rl-1,tk+1)(h)2)+(s(rl+1,tk)-2s(rl,tk)+s(rl-1,tk)(h)2))]+1rl[+(s(rl+1,tk)-s(rl-1,tk)2(h)))12((s(rl+1,tk+1)-s(rl-1,tk+1)2(h))+(s(rl+1,tk)-s(rl-1,tk)2(h)))]+12λ2(s(rl,tk+1)+s(rl,tk)).For simplicity, let us put(12)bjl,k+1=(j+1)1-αlk+1-(j)1-αlk+1,Tlk+1=Γ(2-αlk+1)ταlk+1Sh2T,Glk+1=Γ(2-αlk+1)ταlk+1Sh,λjl,k+1=bj-1l,k+1-bjl,k+1,Γlk+1=Γ(2-αlk+1)ταlk+12λ2.Equation (11) becomes
(13)slk+1(1+2Tlk+1-Γlk+1)=sl+1k+1(Tlk+1+Glk+1rl)+sl-1k+1(Tlk+1-Glk+1rl)+sl+1k(Tlk+1-Glk+1rl)+slk+1(1+2Tlk+1)+∑j=1k(slk+1-j-slk-j)λjl,k+1Glk+1+slk(1-Γlk+1+Tlk+1).


## 4. Stability Analysis of the Crank-Nicolson Scheme

In this section, we will analyze the stability conditions of the Crank-Nicolson scheme for the Hantush equation for a deformable aquifer.LetΥlk=slk-Slk, where Slk is the approximate solution at the point (xl,tk),(k=1,2,…,N,l=1,2,…,M-1), and in addition, Υk=[Υ1k,Υ2k,…,ΥM-1k]T, and the function ζk(x) is chosen to be
(14)Υk(x)={Υlkifxl-h2<x≤xl+h2,l=1,2,…,M-10ifL-h2<x≤L.Then, the functionΥk(x) can be expressed in Fourier series as follows:
(15)Υk(x)=∑m=-∞m=∞δm(m)exp[2iπmkL],δk(x)=1L∫0LΥk(x)exp[2iπmxL]dx.It was established by [15] that
(16)∥Υ2∥22=∑m=-∞m=∞∥δk(m)∥2.Observe that for allk,l≥1,  0≤1-αlk+1<1, and in addition, according to the problem in point, the velocity seepage v, the dispersion coefficient D, the retardation factor R, and the radioactive decay constant λ are positive constants. Then the following properties of the coefficients Tlk+1,Glk+1,λjl,k+1blk+1, and Γlk+1can be established.(1)
Glk+1,Tlk+1,  and  Γlk+1 are positive for all l=1,2,…,M-1,(2)
0<λjl,k≤λj-1l,k≤1 for all l=1,2,…,M-1,(3)
0≤bjl,k≤1,∑j=0k-1bj+1l,k+1=1-λkl,k+1 for all l=1,2,…,M-1.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/rl)→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:
(17)Υlk+1(1+2Tlk+1-Γlk+1)=Υl+1k+1(Tlk+1)+Υl-1k+1(Tlk+1)+Υl+1k(Tlk+1)+Υlk(1+2Tlk+1-Γlk+1)+∑j=1k(Υlk+1-j-Υlk-j)λjl,k+1Glk+1.If we assume thatΥlk in (17) can be put in the delta-exponential form as follows [21]:
(18)Υlk=δkexp[iφlk],
where φ is a real spatial wave number, now replacing the above equation (18) in (17) we obtain
(19)[1+4Tl1sin2(φh2)-2sin2(φh2)Γl1]δ1=[1-4Tl1sin2(φh2)-2sin2(φh2)Γl1]δ0fork=0,[1+4Tl1+ksin2(φh2)-2sin2(φh2)Γl1+k]δk+1=[1-4Tl1+ksin2(φh2)-λ1l,k+1-2sin2(φh2)Γl1+k]δk+∑j=0k-1λj+1l,k+1δk-j+λkl,k+1δ0,fork=1,2,…,N-1.Equation (19) can be written in the following form:
(20)δ1=[1-4Tl1sin2(φh/2)-2sin2(φh/2)Γl1]δ0[1+4Tl1sin2(φh/2)-2sin2(φh/2)Γl1],δk+1=(∑j=0k-1λj+1l,k+1δk-j+λkl,k+1δ0[1-4Tl1+ksin2(φh2)-e1l,k+1-2sin2(φh2)Γl1+k]δk+∑j=0k-1λj+1l,k+1δk-j+λkl,k+1δ0)×[1+4Tl1+ksin2(φh2)-2sin2(φh2)Γl1+k]-1.Our next concern here is to show that for allk=1,2,…,N-1 the solution of (19) satisfies the following condition:
(21)|δk|<|δ0|.To achieve this, we make use of the recurrence technique on the natural numberk.Fork=1, and remembering that dlk+1, blk+1 are positive for all l=1,2,…,M-1, then we obtain
(22)|δ1||δ0|=|[1-4Tl1sin2(φh/2)-2sin2(φh/2)Γl1][1+4Tl1sin2(φh/2)-2sin2(φh/2)Γl1]|<1.Assuming that form=2,3,…,k the property is verified, then
(23)|δk+1|=|(∑j=0k-1λj+1l,k+1δk-j+λkl,k+1[1-4Tl1+ksin2(φh2)-e1l,k+1-2sin2(φh2)Γl1+k]δk+∑j=0k-1λj+1l,k+1δk-j+λkl,k+1δ0)×[1+4Tl1+ksin2(φh2)-2sin2(φh2)Γl1+k]-1(∑j=0k-1λj+1l,k+1δk-j+λkl,k+1[1-4Tl1+ksin2(φh2)-e1l,k+1|.Making use of the triangular inequality, we obtain(24)|δk+1|≤(|∑j=0k-1pj+1l,k+1δk-j||1-4bl1+ksin2(φh2)-e1l,k+1-2sin2(φh2)Γl1+k||δk|+|∑j=0k-1pj+1l,k+1δk-j|+|ekl,k+1δ0|)×|1+4Tl1+ksin2(φh2)-2sin2(φh2)Γl1+k|-1.Using the recurrence hypothesis, we have(25)|δk+1|≤((|∑j=0k-1λj+1l,k+1||1-4Tl1+ksin2(φh2)-2sin2(φh2)Γl1+k|+|∑j=0k-1λj+1l,k+1|)×|1+4Tl1+ksin2(φh2)-2sin2(φh2)Γl1+k|-1(|∑j=0k-1λj+1l,k+1||1-4Tl1+ksin2(φh2)-2sin2(φh2)Γl1+k|)|δ0|,|δk+1|≤(-2sin2(φh2)Γl1+k|-1|1+4Tl1+ksin2(φh2)-2sin2(φh2)Γl1+k|×|1+4Tl1+ksin2(φh2)-2sin2(φh2)Γl1+k|-1)|δ0|,|δk+1|≤|δ0|
which completes the proof that starts at (21) and ends at (25).


## 5. Convergence Analysis of the Crank-Nicolson Scheme

If we assume thats(rl,tk)(l=1,2,…,M,k=1,2,…,N-1) is the exact solution of our problem at the point (rl,tk), by letting Hlk=s(rl,tk)-Slk and Hk=(0,H1k,H2k,…,HM-1k) substituting this in (17), we obtain
(26)Hl1(1+2Tl1-Γl1)-Hl+11(Tl1)-Hl-11(Tl1)=Rl1fork=0,Hl1+k(1+2Tl1+k-Γl1+k)-Hl+1k+1(Tl1+k)-Hl-11+k(Tl1+l)=Rl1+k+∑j=0k-1Hlk-jλj+1l,k+1fork≥1.Here,(27)Rl1+k=s(rl,tk+1)-∑j=0k-1s(rl,tk-j)λj+1l,k+1+b1l,k+1s(rl,t0)+(Γl1+k-Tl1+k)[s(rl+1,tk+1)-2s(rl,tk+1)+s(rl-1,tk+1)].From (8) and (9), we have
(28)∂2s(rl,tk+1)∂r2+h2V1=12((s(rl+1,tk+1)-2s(rl,tk+1)+s(rl-1,tk+1))h2+(s(rl+1,tk)-2s(rl,tk)+s(rl-1,tk))h2),∂αlk+1s(rl,tk+1)∂tαlk+1+τV2=τ-αlk+1Γ(2-αlk+1)(∑j=1k[s(rl,tk+1-j)-s(rl,tk-j)]λjl,ks(rl,tk+1)-s(rl,tk)+∑j=1k[s(rl,tk+1-j)-s(rl,tk-j)]λjl,k).From the above, we have that(29)Rlk+1≤K(τ1+αlk+1+h2ταlk),
where K1,K2, and K 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].Lemma 1.
The following inequality∥Hk+1∥∞≤K(τ1+αlk+1+h2ταlk)(Hjl,k+1)-1 is true for (k=0,1,2,…,N-1)  where ∥Hk∥∞=max1≤l≤M-1(Hk), K is a constant. In addition,
(30)ϕϵ(a,b)={min1≤l≤M-1αlk+1ifτ<1max1≤l≤M-1αlk+1ifτ>1.This can be achieved via the recurrence technique on the natural numberk.Whenk=0, we have the following:
(31)|Hl1|≤(Tl1)|Hl+11|+(Tl1)|Hl-11|=|Fl1|≤V(τ1+αlk+1+h2ταlk)(λjl,k+1)-1.Now suppose that∥Hi+1∥∞≤K(τ1+αli+1+h2ταli)(λjl,i+1)-1,i=1,…,N-2. Then,
(32)|Hl1+k|≤|Hl1+k(1+2Tl1+k-Γl1+k)-Hl+1k+1(Tl1+k)-Hl-11+k(Tl1+l)|≤(Tl1+k)|Hl+1k+1|+(Tl1+k)|Hl-1k+1|+(1+2Tl1+k-Γl1+k)|Hlk+1|=|Rlk+1+∑i=1k(Hlk-i)λjl,k+1|≤|Rlk+1|+∑i=1k|Hlk-i|λjl,k+1≤K(τ1+αlk+1+h2ταlk)+∑i=1k∥Hlk-i∥∞λjl,k+1≤K(τ1+αlk+1+h2ταlk)(λjl,k+1+λ0l,k+1-λjl,k+1)×(λjl,k+1)-1≤V(τ1+αlk+1+h2ταlk)(λ0l,k+1)(λjl,k+1)-1≤V(τ1+αlk+1+h2ταlk)(λjl,k+1)-1,
which completes the proof of Lemma 1.Theorem 2.
The Crank-Nicolson scheme is convergent, and there exists a positive constantV such that
(33)|Slk-s(xl,tk)|≤K(τ+h2),l=1,2,…,M-1,k=1,2,…,N.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.


## 6. 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:(34)SDtα(x,t)s(r,t)=TDrrs(r,t)+1rDrs(r,t)+s(r,t)λ2,α(x,t)=1-sin(x+y),s(r,0)=s0,limr→∞s(r,t)=s0,Q=2πn/2Γ(n/2)rbn-1Kd3-n∂rs(rb,t).
Figure 2 shows the numerical simulation for an aquifer thickness of 1000 feet, a transmissivity T=10000 feet/day, a hydraulic conductivity of the main aquifer of K=8 feet/day, and a hydraulic conductivity of the leaky Kl=0.4 feet/day. A storativity S=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 r=400 feet. The black shows the simulation for a distance of r=200 feet, and finally the blue is the simulation for r=3650 feet. Figure 3 shows the numerical simulation of the water flowing in the variable leaky aquifer.Figure 2
Numerical comparison of the observed data and solution of the modified equation for different values of time.Figure 3
Numerical simulation of the water flowing in the leaky aquifer.


## 7. 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.



---

*Source: 691060-2013-05-27.xml*