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

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.

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

Let 𝑓∶ℝ→ℝ,𝑥→𝑓(𝑥) denote a continuous but necessary differentiable; let 𝛼(𝑥) be a continuous function in (0, 1]. Then its variational order differential is defined as 𝐷₀^𝛼(𝑥)𝑓(𝑥)=1Γ(1−𝛼(𝑥))^𝑥₀(𝑥−𝑡)^{−𝛼(𝑡)}𝑑𝑓(𝑡)𝑑𝑡𝑑𝑡. (1)

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 𝑆𝐷_𝑡Φ(𝑟,𝑡)=𝑇𝐷_𝑟𝑟Φ()+1𝑟,𝑡𝑟𝐷_𝑟Φ(𝑟,𝑡). (2)

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: 𝑆𝐷_𝑡𝑠(𝑟,𝑡)=𝑇𝐷_𝑟𝑟𝑠()+1𝑟,𝑡𝑟𝐷_𝑟𝑠()+𝑠()𝑟,𝑡𝑟,𝑡𝜆², (3) where 𝑠 is the drawdown or change in the level of water; 𝑆 is the specific storativity of the aquifer, and 𝑇 is the transmissivity: 𝜆²=𝐵𝐾^′, (4) 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: 𝑠(𝑟,0)=Φ₀,lim_𝑟→∞Φ(𝑟,𝑡)=Φ₀,𝑄=2𝜋^𝑛/2Γ()𝑟𝑛/2^𝑏𝑛−1𝐾𝑑^3−𝑛𝜕_𝑟𝑠𝑟_𝑏,𝑡. (5)

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 𝑆𝐷_𝑡^{𝛼(𝑥,𝑡)}𝑠(𝑟,𝑡)=𝑇𝐷_𝑟𝑟𝑠()+1𝑟,𝑡𝑟𝐷_𝑟𝑠()+𝑠()𝑟,𝑡𝑟,𝑡𝜆²,0<𝛼(𝑥,𝑡)≤1. (6)

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 𝑥_𝑙=𝑙ℎ,  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 𝜕𝑠=1𝜕𝑟2𝑠𝑟_𝑙+1,𝑡_𝑘+1−Φ𝑟_𝑙−1,𝑡_𝑘+12(ℎ)+𝑠𝑟_𝑙+1,𝑡_𝑘−𝑠𝑟_𝑙−1,𝑡_𝑘2(ℎ)(ℎ),𝜕+𝑂₂𝑠𝜕𝑟²=12𝑠𝑟_𝑙+1,𝑡_𝑘+1−2𝑠𝑟_𝑙,𝑡_𝑘+1+𝑠𝑟_𝑙−1,𝑡_𝑘+1(ℎ)²+𝑠𝑟_𝑙+1,𝑡_𝑘−2𝑠𝑟_𝑙,𝑡_𝑘+𝑠𝑟_𝑙−1,𝑡_𝑘(ℎ)²+𝑂ℎ²,1𝑠=2𝑠𝑟_𝑙,𝑡_𝑘+1+𝑠𝑟_𝑙,𝑡_𝑘. (7)(8)(9)

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}. (10)

Now replacing (7), (8), (9), and (10) in (6), we obtain the following: 𝑆𝜏^{−𝛼𝑙𝑘+1}Γ2−𝛼_𝑙^𝑘+1𝑠𝑟_𝑙,𝑡_𝑘+1−𝑠𝑟_𝑙,𝑡_𝑘+^𝑘_𝑗=1𝑠𝑟_𝑙,𝑡_{𝑘+1−𝑗}−𝑠𝑟_𝑙,𝑡_𝑘−𝑗×𝑗+1^{1−𝛼𝑙𝑘+1}−𝑗^{1−𝛼𝑙𝑘+1}1=𝑇2𝑠𝑟_𝑙+1,𝑡_𝑘+1−2𝑠𝑟_𝑙,𝑡_𝑘+1+𝑠𝑟_𝑙−1,𝑡_𝑘+1(ℎ)²+𝑠𝑟_𝑙+1,𝑡_𝑘−2𝑠𝑟_𝑙,𝑡_𝑘+𝑠𝑟_𝑙−1,𝑡_𝑘(ℎ)²+1𝑟_𝑙12𝑠𝑟_𝑙+1,𝑡_𝑘+1−𝑠𝑟_𝑙−1,𝑡_𝑘+12(ℎ)+𝑠𝑟_𝑙+1,𝑡_𝑘−𝑠𝑟_𝑙−1,𝑡_𝑘2(ℎ)+12𝜆²𝑠𝑟_𝑙,𝑡_𝑘+1+𝑠𝑟_𝑙,𝑡_𝑘. (11)

For simplicity, let us put 𝑏_𝑗^𝑙,𝑘+1=𝑗+1^{1−𝛼𝑙𝑘+1}−𝑗^{1−𝛼𝑙𝑘+1},𝑇_𝑙^𝑘+1=Γ2−𝛼_𝑙^𝑘+1𝜏^{𝛼𝑙𝑘+1}𝑆ℎ²𝑇,𝐺_𝑙^𝑘+1=Γ2−𝛼_𝑙^𝑘+1𝜏^{𝛼𝑙𝑘+1},𝜆𝑆ℎ_{𝑗𝑙,𝑘+1}=𝑏^𝑙,𝑘+1_𝑗−1−𝑏_𝑗^𝑙,𝑘+1,Γ_𝑙^𝑘+1=Γ2−𝛼_𝑙^𝑘+1𝜏^{𝛼𝑙𝑘+1}2𝜆². (12)

Equation (11) becomes 𝑠_𝑙^𝑘+11+2𝑇_𝑙^𝑘+1−Γ_𝑙^𝑘+1=𝑠^𝑘+1_𝑙+1𝑇_𝑙^𝑘+1+𝐺_𝑙^𝑘+1𝑟_𝑙+𝑠^𝑘+1_𝑙−1𝑇_𝑙^𝑘+1−𝐺_𝑙^𝑘+1𝑟_𝑙+𝑠^𝑘_𝑙+1𝑇_𝑙^𝑘+1−𝐺_𝑙^𝑘+1𝑟_𝑙+𝑠_𝑙^𝑘+11+2𝑇_𝑙^𝑘+1+^𝑘_𝑗=1𝑠_𝑙^{𝑘+1−𝑗}−𝑠_𝑙^𝑘−𝑗𝜆_𝑗^𝑙,𝑘+1𝐺_𝑙^𝑘+1+𝑠^𝑘_𝑙1−Γ_𝑙^𝑘+1+𝑇_𝑙^𝑘+1. (13)

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 Υ^𝑘_𝑙=𝑠^𝑘_𝑙−𝑆^𝑘_𝑙, where 𝑆^𝑘_𝑙 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<𝑥≤𝐿. (14)

Then, the function Υ^𝑘(𝑥) can be expressed in Fourier series as follows: Υ^𝑘(𝑥)=^𝑚=∞_𝑚=−∞𝛿_𝑚(𝑚)exp2𝑖𝜋𝑚𝑘𝐿,𝛿_𝑘(𝑥)=1𝐿^𝐿₀Υ^𝑘(𝑥)exp2𝑖𝜋𝑚𝑥𝐿𝑑𝑥. (15)

It was established by [15] that Υ²²₂=^𝑚=∞_𝑚=−∞𝛿_𝑘(𝑚)². (16)

Observe that for all 𝑘,𝑙≥1,  0≤1−𝛼_𝑙^𝑘+1<1, and 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𝑏_𝑙^𝑘+1, and Γ_𝑙^𝑘+1can be established.(1)𝐺_𝑙^𝑘+1,𝑇_𝑙^𝑘+1,  and  Γ_𝑙^𝑘+1 are positive for all 𝑙=1,2,…,𝑀−1,(2)0<𝜆_𝑗^𝑙,𝑘≤𝜆^𝑙,𝑘_𝑗−1≤1 for all 𝑙=1,2,…,𝑀−1,(3)0≤𝑏_𝑗^𝑙,𝑘≤1,∑^𝑘−1_𝑗=0𝑏^𝑙,𝑘+1_𝑗+1=1−𝜆_𝑘^𝑙,𝑘+1 for all 𝑙=1,2,…,𝑀−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/𝑟_𝑙)→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: Υ_𝑙^𝑘+11+2𝑇_𝑙^𝑘+1−Γ_𝑙^𝑘+1=Υ^𝑘+1_𝑙+1𝑇_𝑙^𝑘+1+Υ^𝑘+1_𝑙−1𝑇_𝑙^𝑘+1+Υ^𝑘_𝑙+1𝑇_𝑙^𝑘+1+Υ^𝑘_𝑙1+2𝑇_𝑙^𝑘+1−Γ_𝑙^𝑘+1+^𝑘_𝑗=1Υ_𝑙^{𝑘+1−𝑗}−Υ_𝑙^𝑘−𝑗𝜆_𝑗^𝑙,𝑘+1𝐺_𝑙^𝑘+1. (17)

If we assume that Υ^𝑘_𝑙 in (17) can be put in the delta-exponential form as follows [21]: Υ^𝑘_𝑙=𝛿_𝑘exp𝑖𝜑𝑙𝑘, (18) where 𝜑 is a real spatial wave number, now replacing the above equation (18) in (17) we obtain 1+4𝑇¹_𝑙sin²𝜑ℎ2−2sin²𝜑ℎ2Γ¹_𝑙𝛿₁=1−4𝑇¹_𝑙sin²𝜑ℎ2−2sin²𝜑ℎ2Γ¹_𝑙𝛿₀for𝑘=0,1+4𝑇_𝑙^1+𝑘sin²𝜑ℎ2−2sin²𝜑ℎ2Γ_𝑙^1+𝑘𝛿_𝑘+1=1−4𝑇_𝑙^1+𝑘sin²𝜑ℎ2−𝜆₁^𝑙,𝑘+1−2sin²𝜑ℎ2Γ_𝑙^1+𝑘𝛿_𝑘+^𝑘−1_𝑗=0𝜆^𝑙,𝑘+1_𝑗+1𝛿_𝑘−𝑗+𝜆_𝑘^𝑙,𝑘+1𝛿₀,for𝑘=1,2,…,𝑁−1. (19)

Equation (19) can be written in the following form: 𝛿₁=1−4𝑇¹_𝑙sin²𝜑ℎ/2−2sin²𝜑ℎ/2Γ¹_𝑙𝛿₀1+4𝑇¹_𝑙sin²𝜑ℎ/2−2sin²𝜑ℎ/2Γ¹_𝑙,𝛿_𝑘+1=1−4𝑇_𝑙^1+𝑘sin²𝜑ℎ2−𝑒₁^𝑙,𝑘+1−2sin²𝜑ℎ2Γ_𝑙^1+𝑘𝛿_𝑘+^𝑘−1_𝑗=0𝜆^𝑙,𝑘+1_𝑗+1𝛿_𝑘−𝑗+𝜆_𝑘^𝑙,𝑘+1𝛿₀×1+4𝑇_𝑙^1+𝑘sin²𝜑ℎ2−2sin²𝜑ℎ2Γ_𝑙^1+𝑘−1. (20)

Our next concern here is to show that for all 𝑘=1,2,…,𝑁−1 the solution of (19) satisfies the following condition: 𝛿_𝑘<𝛿₀. (21)

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−2sin²𝜑ℎ/2Γ¹_𝑙1+4𝑇¹_𝑙sin²𝜑ℎ/2−2sin²𝜑ℎ/2Γ¹_𝑙<1. (22)

Assuming that for 𝑚=2,3,…,𝑘 the property is verified, then 𝛿_𝑘+1=1−4𝑇_𝑙^1+𝑘sin²𝜑ℎ2−𝑒₁^𝑙,𝑘+1−2sin²𝜑ℎ2Γ_𝑙^1+𝑘𝛿_𝑘+^𝑘−1_𝑗=0𝜆^𝑙,𝑘+1_𝑗+1𝛿_𝑘−𝑗+𝜆_𝑘^𝑙,𝑘+1𝛿₀×1+4𝑇_𝑙^1+𝑘sin²𝜑ℎ2−2sin²𝜑ℎ2Γ_𝑙^1+𝑘−1. (23)

Making use of the triangular inequality, we obtain 𝛿_𝑘+1≤1−4𝑏_𝑙^1+𝑘sin²𝜑ℎ2−𝑒₁^𝑙,𝑘+1−2sin²𝜑ℎ2Γ_𝑙^1+𝑘𝛿_𝑘+^𝑘−1_𝑗=0𝑝^𝑙,𝑘+1_𝑗+1𝛿_𝑘−𝑗+𝑒_𝑘^𝑙,𝑘+1𝛿₀×1+4𝑇_𝑙^1+𝑘sin²𝜑ℎ2−2sin²𝜑ℎ2Γ_𝑙^1+𝑘−1. (24)

Using the recurrence hypothesis, we have 𝛿_𝑘+1≤1−4𝑇_𝑙^1+𝑘sin²𝜑ℎ2−2sin²𝜑ℎ2Γ_𝑙^1+𝑘+^𝑘−1_𝑗=0𝜆^𝑙,𝑘+1_𝑗+1×1+4𝑇_𝑙^1+𝑘sin²𝜑ℎ2−2sin²𝜑ℎ2Γ_𝑙^1+𝑘−1𝛿₀,𝛿_𝑘+1≤1+4𝑇_𝑙^1+𝑘sin²𝜑ℎ2−2sin²𝜑ℎ2Γ_𝑙^1+𝑘×1+4𝑇_𝑙^1+𝑘sin²𝜑ℎ2−2sin²𝜑ℎ2Γ_𝑙^1+𝑘−1𝛿₀,𝛿_𝑘+1≤𝛿₀ (25) which completes the proof that starts at (21) and ends at (25).

5. 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) substituting this in (17), we obtain 𝐻¹_𝑙1+2𝑇¹_𝑙−Γ¹_𝑙−𝐻¹_𝑙+1𝑇¹_𝑙−𝐻¹_𝑙−1𝑇¹_𝑙=𝑅¹_𝑙𝐻for𝑘=0,_𝑙1+𝑘1+2𝑇_𝑙^1+𝑘−Γ_𝑙^1+𝑘−𝐻^𝑘+1_𝑙+1𝑇_𝑙^1+𝑘−𝐻^1+𝑘_𝑙−1𝑇_𝑙^1+𝑙=𝑅_𝑙^1+𝑘+^𝑘−1_𝑗=0𝐻_𝑙^𝑘−𝑗𝜆^𝑙,𝑘+1_𝑗+1for𝑘≥1. (26)

Here, 𝑅_𝑙^1+𝑘=𝑠𝑟_𝑙,𝑡_𝑘+1−^𝑘−1_𝑗=0𝑠𝑟_𝑙,𝑡_𝑘−𝑗𝜆^𝑙,𝑘+1_𝑗+1+𝑏₁^𝑙,𝑘+1𝑠𝑟_𝑙,𝑡₀+Γ_𝑙^1+𝑘−𝑇_𝑙^1+𝑘𝑠𝑟_𝑙+1,𝑡_𝑘+1−2𝑠𝑟_𝑙,𝑡_𝑘+1+𝑠𝑟_𝑙−1,𝑡_𝑘+1. (27)

From (8) and (9), we have 𝜕²𝑠𝑟_𝑙,𝑡_𝑘+1𝜕𝑟²+ℎ²𝑉₁=12𝑠𝑟_𝑙+1,𝑡_𝑘+1−2𝑠𝑟_𝑙,𝑡_𝑘+1+𝑠𝑟_𝑙−1,𝑡_𝑘+1ℎ²+𝑠𝑟_𝑙+1,𝑡_𝑘−2𝑠𝑟_𝑙,𝑡_𝑘+𝑠𝑟_𝑙−1,𝑡_𝑘ℎ²,𝜕^{𝛼𝑙𝑘+1}𝑠𝑟_𝑙,𝑡_𝑘+1𝜕𝑡^{𝛼𝑙𝑘+1}+𝜏𝑉₂=𝜏^{−𝛼𝑙𝑘+1}Γ2−𝛼_𝑙^𝑘+1𝑠𝑟_𝑙,𝑡_𝑘+1−𝑠𝑟_𝑙,𝑡_𝑘+^𝑘_𝑗=1𝑠𝑟_𝑙,𝑡_{𝑘+1−𝑗}−𝑠𝑟_𝑙,𝑡_𝑘−𝑗𝜆_𝑗^𝑙,𝑘. (28)

From the above, we have that 𝑅_𝑙^𝑘+1≤𝐾𝜏^{1+𝛼𝑙𝑘+1}+ℎ²𝜏^𝛼𝑘𝑙, (29) where 𝐾₁, 𝐾₂, 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].

Lemma 1. The following inequality ‖𝐻^𝑘+1‖_∞≤𝐾(𝜏^{1+𝛼𝑙𝑘+1}+ℎ²𝜏^𝛼𝑘𝑙)(𝐻_𝑗^𝑙,𝑘+1)⁻¹ is true for (𝑘=0,1,2,…,𝑁−1) where ‖𝐻^𝑘‖_∞=max_{1≤𝑙≤𝑀−1}(𝐻^𝑘), 𝐾 is a constant. In addition, 𝜙^𝜖(𝑎,𝑏)=min_{1≤𝑙≤𝑀−1}𝛼_𝑙^𝑘+1if𝜏<1max_{1≤𝑙≤𝑀−1}𝛼_𝑙^𝑘+1if𝜏>1. (30)

This can be achieved via the recurrence technique on the natural number 𝑘.

When 𝑘=0, we have the following: 𝐻¹_𝑙≤𝑇¹_𝑙𝐻¹_𝑙+1+𝑇¹_𝑙𝐻¹_𝑙−1=𝐹¹_𝑙≤𝑉𝜏^{1+𝛼𝑙𝑘+1}+ℎ²𝜏^𝛼𝑘𝑙𝜆_𝑗^{𝑙,𝑘+1−1}. (31)

Now suppose that ‖𝐻^𝑖+1‖_∞≤𝐾(𝜏^{1+𝛼𝑙𝑖+1}+ℎ²𝜏^𝛼𝑖𝑙)(𝜆_𝑗^𝑙,𝑖+1)⁻¹,𝑖=1,…,𝑁−2. Then, 𝐻_𝑙^1+𝑘≤𝐻_𝑙^1+𝑘1+2𝑇_𝑙^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}, (32) which completes the proof of Lemma 1.

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

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: 𝑆𝐷_𝑡^{𝛼(𝑥,𝑡)}𝑠(𝑟,𝑡)=𝑇𝐷_𝑟𝑟𝑠()+1𝑟,𝑡𝑟𝐷_𝑟𝑠()+𝑠()𝑟,𝑡𝑟,𝑡𝜆²,𝛼(),𝑠()𝑥,𝑡=1−sin𝑥+𝑦𝑟,0=𝑠₀,lim_𝑟→∞𝑠(𝑟,𝑡)=𝑠₀,𝑄=2𝜋^𝑛/2Γ()𝑟𝑛/2^𝑏𝑛−1𝐾𝑑^3−𝑛𝜕_𝑟𝑠𝑟_𝑏,𝑡. (34) Figure 2 shows the numerical simulation for an aquifer thickness of 1000 feet, a transmissivity 𝑇=10000 feet/day, a 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.

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.