Abstract and Applied Analysis

Volume 2013 (2013), Article ID 691060, 8 pages

http://dx.doi.org/10.1155/2013/691060

Research Article

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

Abdon Atangana¹ and S. C. Oukouomi Noutchie²

¹Institute for Groundwater Studies, University of the Free State, P.O. Box 399, Bloemfontein, South Africa

²Department of Mathematical Sciences, North-West University, Mafikeng Campus, Mmabatho 2735, South Africa

Received 21 March 2013; Revised 30 April 2013; Accepted 7 May 2013

Academic Editor: Bashir Ahmad

Copyright © 2013 Abdon Atangana and S. C. Oukouomi Noutchie. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

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

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

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 ,    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

The Crank-Nicolson scheme for the VO time-fractional diffusion model can be stated as follows:

Now replacing (7), (8), (9), and (10) in (6), we obtain the following:

For simplicity, let us put

Equation (11) becomes

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 , and in addition, , and the function is chosen to be

Then, the function can be expressed in Fourier series as follows:

It was established by [15] that

Observe that for all ,  , 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 , and can be established.

(1),  and   are positive for all ,(2) for all ,(3) for all .

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   [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 delta-exponential 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 the solution of (19) satisfies the following condition:

To achieve this, we make use of the recurrence technique on the natural number .

For , and remembering that , are positive for all , then we obtain

Assuming that for the property is verified, then

Making use of the triangular inequality, we obtain

Using the recurrence hypothesis, we have

which completes the proof that starts at (21) and ends at (25).

5. Convergence Analysis of the Crank-Nicolson Scheme

If we assume that is the exact solution of our problem at the point , by letting and substituting this in (17), we obtain

Here,

From (8) and (9), we have

From the above, we have that

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 is true for where , is a constant. In addition,

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

When , we have the following:

Now suppose that . Then,

which completes the proof of Lemma 1.

Theorem 2. The Crank-Nicolson scheme is convergent, and there exists a positive constant such that

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:

Figure 2 shows the numerical simulation for an aquifer thickness of 1000 feet, a transmissivity  feet/day, a hydraulic conductivity of the main aquifer of  feet/day, and a hydraulic conductivity of the leaky  feet/day. A storativity , the leaky factor is 0.00081/day, and finally the flow rate  feet/day. The red shows the simulation for a distance of feet. The black shows the simulation for a distance of feet, and finally the blue is the simulation for 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.

References

G. P. Kruseman and N. A. de Ridder, Analysis and Evaluation of Pumping Test Data, International Institute for Land Reclamation and Improvement, Wageningen, The Netherlands, 2nd edition, 1990.

C. V. Theis, “The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using ground-water storage,” Transactions, American Geophysical Union, vol. 16, no. 2, pp. 519–524, 1935.

M. S. Hantush, “Analysis of data from pumping tests in leaky aquifers,” Transactions, American Geophysical Union, vol. 37, no. 6, pp. 702–714, 1956.

M. S. Hantush and C. E. Jacob, “Non-steady radial flow in an infinite leaky aquifer,” Transactions, American Geophysical Union, vol. 36, no. 1, pp. 95–100, 1955.

S. Umarov and S. Steinberg, “Variable order differential equations with piecewise constant order-function and diffusion with changing modes,” Journal of Analysis and Its Applications, vol. 28, no. 4, pp. 431–450, 2009.

B. Ross and S. Samko, “Fractional integration operator of variable order in the Hölder spaces ,” International Journal of Mathematics and Mathematical Sciences, vol. 18, no. 4, pp. 777–788, 1995.

A. Atangana and J. F. Botha, “Generalized groundwater flow equation using the concept of variable order derivative,” Boundary Value Problems, vol. 2013, article 53, 2013.

A. Atangana and A. Kılıçman, “A possible generalization of acoustic wave equation using the concept of perturbed derivative order,” Mathematical Problems in Engineering, vol. 2013, Article ID 696597, 6 pages, 2013.

Y. Zhang, “A finite difference method for fractional partial differential equation,” Applied Mathematics and Computation, vol. 215, no. 2, pp. 524–529, 2009.

C. Tadjeran, M. M. Meerschaert, and H.-P. Scheffler, “A second-order accurate numerical approximation for the fractional diffusion equation,” Journal of Computational Physics, vol. 213, no. 1, pp. 205–213, 2006.

M. M. Meerschaert and C. Tadjeran, “Finite difference approximations for fractional advection-dispersion flow equations,” Journal of Computational and Applied Mathematics, vol. 172, no. 1, pp. 65–77, 2004.

W. H. Deng, “Numerical algorithm for the time fractional Fokker-Planck equation,” Journal of Computational Physics, vol. 227, no. 2, pp. 1510–1522, 2007.

C. P. Li, A. Chen, and J. J. Ye, “Numerical approaches to fractional calculus and fractional ordinary differential equation,” Journal of Computational Physics, vol. 230, no. 9, pp. 3352–3368, 2011.

C.-M. Chen, F. Liu, I. Turner, and V. Anh, “A Fourier method for the fractional diffusion equation describing sub-diffusion,” Journal of Computational Physics, vol. 227, no. 2, pp. 886–897, 2007.

Y. Q. Chen and K. L. Moore, “Discretization schemes for fractional-order differentiators and integrators,” IEEE Transactions on Circuits and Systems I, vol. 49, no. 3, pp. 363–367, 2002.

S. B. Yuste and L. Acedo, “An explicit finite difference method and a new von Neumann-type stability analysis for fractional diffusion equations,” SIAM Journal on Numerical Analysis, vol. 42, no. 5, pp. 1862–1874, 2005.

I. Podlubny, A. Chechkin, T. Skovranek, Y. Q. Chen, and B. M. Vinagre Jara, “Matrix approach to discrete fractional calculus. II. Partial fractional differential equations,” Journal of Computational Physics, vol. 228, no. 8, pp. 3137–3153, 2009.

E. Hanert, “On the numerical solution of space-time fractional diffusion models,” Computers and Fluids, vol. 46, pp. 33–39, 2011.

Y. Zhang, “A finite difference method for fractional partial differential equation,” Applied Mathematics and Computation, vol. 215, no. 2, pp. 524–529, 2009.

R. Lin, F. Liu, V. Anh, and I. Turner, “Stability and convergence of a new explicit finite-difference approximation for the variable-order nonlinear fractional diffusion equation,” Applied Mathematics and Computation, vol. 212, no. 2, pp. 435–445, 2009.

H. Sun, W. Chen, C. Li, and Y. Chen, “Finite difference schemes for variable-order time fractional diffusion equation,” International Journal of Bifurcation and Chaos in Applied Sciences and Engineering, vol. 22, no. 4, Article ID 1250085, 16 pages, 2012.

J. Crank and P. Nicolson, “A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type,” Mathematical Proceedings of the Cambridge Philosophical Society, vol. 43, no. 1, pp. 50–67, 1947.

C. P. Li and C. X. Tao, “On the fractional Adams method,” Computers and Mathematics with Applications, vol. 58, no. 8, pp. 1573–1588, 2009.