SpringerOpen Newsletter

Receive periodic news and updates relating to SpringerOpen.

This article is part of the series Proceedings of the International Congress in Honour of Professor Hari M. Srivastava.

Open Access Research

A generalized groundwater flow equation using the concept of variable-order derivative

Abdon Atangana* and Joseph Francois Botha

Author Affiliations

Institute for Groundwater Studies, Faculty of Natural and Agricultural Sciences, University of the Free State, P.O. Box 9301, Bloemfontein, South Africa

For all author emails, please log on.

Boundary Value Problems 2013, 2013:53  doi:10.1186/1687-2770-2013-53

The electronic version of this article is the complete one and can be found online at: http://www.boundaryvalueproblems.com/content/2013/1/53


Received:30 January 2013
Accepted:23 February 2013
Published:14 March 2013

© 2013 Atangana and Botha; licensee Springer

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

In this paper, the groundwater flow equation is generalized using the concept of the variational order derivative. We present a numerical solution of the modified groundwater flow equation with the variational order derivative. We solve the generalized equation with the Crank-Nicholson technique. Numerical methods typically yield approximate solutions to the governing equation through the discretization of space and time and can relax the rigid idealized conditions of analytical models or lumped-parameter models. They can therefore be more realistic and flexible for simulating field conditions. Within the discredited problem domain, the variable internal properties, boundaries, and stresses of the system are approximated. We perform the stability and convergence analysis of the Crank-Nicholson method and complete the paper with some illustrative computational examples and their simulations.

Keywords:
groundwater flow equation; variable order derivative; Crank-Nicholson scheme; stability; convergence

1 Introduction

A problem that arises naturally in groundwater investigations is to choose an appropriate geometry for the geological system in which the flow occurs. This geological formation, through which the groundwater flows, changes in time and space. Figure 1 shows the groundwater flow direction toward the pumping well during the extraction of water from the aquifer. The main effort is to establish a suitable mathematical relation between different observables, or what [1] calls a mathematical model for the observables. An ideal mathematical model should not merely provide links between different observables, but also lead to a better understanding of the phenomenon. The attractiveness of this approach is that the mathematical model could, in principle, be used to investigate the future behavior of a given phenomenon under various conditions. This leads us to the groundwater flow equation. This model is based on the conventional, saturated groundwater flow equation for density-independent flow:

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M1">View MathML</a>

(1.1)

thumbnailFigure 1. Groundwater flow in the borehole during pumping test.

Here <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M2">View MathML</a> is the specific storativity, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M3">View MathML</a> is the hydraulic conductivity tensor of the aquifer, Φ is the piezometric head, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M4">View MathML</a> is the strength of any sources or sinks. As a review of the derivation of Eq. (1.1), we will show the Darcy law

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M5">View MathML</a>

(1.2)

is used as a keystone in the derivation of Eq. (1.1). This law, proposed by Darcy early in the nineteenth century, relies on experimental results obtained from the flow of water through a one-dimensional sand column. Alternatively, Darcy’s law states that the rate of flow through a porous medium is proportional to the loss of head and inversely proportional to the length of the flow path. Note that the specific discharge <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M6">View MathML</a> has the dimensions of velocity. Recent investigations [2] suggest that the flow is also influenced by the geometry of the bedding parallel fractures, the feature that equation (1.1) cannot account for. It is therefore possible that equation (1.1) may not be applicable to the flow in these fractured aquifers. In an attempt to circumvent this problem, Barker [3] introduced the model in which the geometry of the aquifer is regarded as a fractal. Although this model has been applied with reasonable success in the analysis of hydraulic tests from boreholes in Karoo aquifers [4], it introduces parameters for which no sound definition exists in the case of non-integer flow dimensions. Recently [5,6], the concept of a fractional-order derivative was used to generalize the groundwater flow equation. However, it has been found that the constant-order fractional diffusion equations are not capable of characterizing some complex diffusion processes, for instance, diffusion process in aninhomogeneous or heterogeneous medium [7].

In addition, when we consider the diffusion process in a porous medium, if the medium structure or external field changes with time, in this situation, the constant-order fractional diffusion equation model cannot be used to well characterize such a phenomenon [8-19]. This is the case of the groundwater flow problem, the medium through which the flow occurs is heterogeneous and changes with time. Still in some biology diffusion processes, the concentration of particles will determine the diffusion pattern [10,11]. To solve the above problems, the variable-order (VO) fractional diffusion equation models have been suggested for use [12]. This present work is therefore devoted to the discussion underpinning the description of the groundwater flow equation with the variable-order derivative.

2 Modified groundwater flow equation

For the readers that are not acquainted with the concept of the variational order derivative, we start this section by presenting the basic definition of this derivative.

2.1 Variational order differential operator

Let <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M7">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M8">View MathML</a> denote a continuous and necessary differentiable, let <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M9">View MathML</a> be a continuous function in <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M10">View MathML</a>. Then its variational order differential in <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M11">View MathML</a> is defined as

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M12">View MathML</a>

(2.1)

The above derivative is called the Caputo variational order differential operator; in addition, 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, 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 the 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. Theis (1935) [13] 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

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M13">View MathML</a>

(2.2)

The above equation is classified under a parabolic equation. 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

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M14">View MathML</a>

(2.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 the analysis of analytical models. The discussion presented in this section will therefore be devoted to the derivation of a numerical solution to groundwater flow equation (2.3).

Numerical methods yield approximate solutions to the governing equation through the discretization of space and time. 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 idealized conditions of analytical models or lumped-parameter models, and they can therefore be more realistic and flexible for simulating field conditions. The finite difference schemes for constant-order time or space fractional diffusion equations have been widely studied [14-19]. For constant-order time fractional diffusion equations, Chen et al. proposed an implicit difference approximation scheme [20]. Yuste et al. introduced weighted average finite difference methods [21]. Podlubny et al. proposed the matrix approach for fractional diffusion equations [22], and Hanert proposed a flexible numerical scheme for the discretization of the space-time fractional diffusion equation [23]. Recently, Zhuang et al. considered the numerical schemes for VO space fractional advection-dispersion equation [16], Lin et al. investigated the explicit scheme for VO nonlinear space fractional diffusion equation [24]. Before applying the numerical methods, we assume Eq. (2.3) has a unique and sufficiently smooth solution. To establish the numerical schemes for the above equation, we let <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M15">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M16">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M17">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M18">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M19">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M20">View MathML</a>, h is the step and τ is the time size, M and N are grid points.

3.1 Crank-Nicholson scheme [24]

We introduce the Crank-Nicholson scheme as follows. Firstly, the discretization of first- and second-order space derivative is stated as follows:

(3.1)

(3.2)

(3.3)

The Crank-Nicholson scheme for the VO time fractional groundwater flow model can be stated as follows:

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M24">View MathML</a>

(3.4)

Now, replacing equations (3.1), (3.2), (3.3) and (3.4) in (2.3), we obtain the following:

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M25">View MathML</a>

For simplicity, let us put:

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M26">View MathML</a>

(3.5)

Equation (3.5) becomes

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M27">View MathML</a>

(3.6)

4 Stability analysis of the Crank-Nicholson scheme

In this section, we will analyze the stability conditions of the Crank-Nicholson scheme for the generalized groundwater flow equation.

Let <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M28">View MathML</a>. Here <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M29">View MathML</a> is the approximate solution at the point <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M30">View MathML</a> (<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M31">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M32">View MathML</a>) and, in addition, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M33">View MathML</a> and the function <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M34">View MathML</a> is chosen to be

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M35">View MathML</a>

(4.1)

Then the function <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M34">View MathML</a> can be expressed in Fourier series as follows:

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M37">View MathML</a>

(4.2)

It was established by Chen et al.[19] that

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M38">View MathML</a>

(4.3)

Observe that for all <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M39">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M40">View MathML</a>, in addition, according to the problem in point, the storativity S, and transmissivity T are positive constants. Then the following properties of the coefficients <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M41">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M42">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M43">View MathML</a>, and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M44">View MathML</a> can be established:

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M45">View MathML</a>

(4.4)

It is customary in groundwater investigations to choose a point on the centerline 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. In such situations where the distribution of the piezometric head in the aquifer is a decreasing function of the distance from the borehole, the expression <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M46">View MathML</a>. Under this situation, the error committed while approximating the solution of the generalized groundwater flow equation with the Crank-Nicholson scheme can be presented as follows:

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M47">View MathML</a>

(4.5a)

If we assume that <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M48">View MathML</a> in equation (4.1) can be put in the delta-exponential form as follows:

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M49">View MathML</a>

(4.5b)

where ϕ is a real spatial wave number, new replacing the above equation (4.5b) into (4.5a), we obtain

(4.6)

Equation (4.6) can be written in the following form:

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M51">View MathML</a>

(4.7)

Our next concern here is to show that for all <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M52">View MathML</a>, the solution of equation (4.7) satisfies the following condition:

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M53">View MathML</a>

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

For <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M54">View MathML</a> and remembering that <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M55">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M44">View MathML</a> are positive for all <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M57">View MathML</a>, we obtain

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M58">View MathML</a>

(4.8)

Assuming that for <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M59">View MathML</a> the property is verified, we get

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M60">View MathML</a>

(4.9)

Making use of the triangular inequality we obtain

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M61">View MathML</a>

(4.10)

Using the recurrence hypothesis, we have

(4.11)

which this completes the proof.

5 Convergence analysis of the Crank-Nicholson scheme

If we assume that <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M63">View MathML</a> (<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M64','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M64">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M52">View MathML</a>) is the exact solution of our problem at the point <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M66','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M66">View MathML</a>, by letting <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M67','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M67">View MathML</a> and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M68">View MathML</a>, substituting this in equation (4.7), we obtain

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M69','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M69">View MathML</a>

(5.1)

Here,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M70">View MathML</a>

(5.2)

From equations (3.1) and (3.4), we have

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M71">View MathML</a>

From the above, we have that

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M72','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M72">View MathML</a>

(5.3)

where <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M73','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M73">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M74','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M74">View MathML</a>, and K are constants. Taking into account the Caputo-type fractional derivative, the detailed error analysis on the above schemes can refer to the work by Diethelm et al.[25] and further work by Li and Tao [26].

Lemma 1<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M75','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M75">View MathML</a>is true for (<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M76','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M76">View MathML</a>), where<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M77','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M77">View MathML</a>, Kis a constant. In addition,

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M78','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M78">View MathML</a>

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

When <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M79','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M79">View MathML</a>, we have the following:

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M80','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M80">View MathML</a>

(5.4)

Now suppose that <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M81','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M81">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M82','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M82">View MathML</a>. Then

(5.5)

which completes the proof.

Theorem 1The Crank-Nicholson scheme is convergent, and there exists a positive constantVsuch that

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M84','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M84">View MathML</a>

(5.6)

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

6 Numerical results

An image is worth ten thousand words; therefore, we devote this section to the numerical simulations of the solution of the generalized groundwater flow equation. The parameters used in the simulation are given as <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M85','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M85">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M86','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M86">View MathML</a>, <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M87','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M87">View MathML</a>, and <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M88','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M88">View MathML</a>.

Example

<a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M89','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M89">View MathML</a>

(6.1)

The numerical simulations of the approximate solution are presented in figures. First, Figure 2 shows the solution as a function of time and space. Figure 3 shows the solution as a function of time for a fixed distance. Figure 4 shows the solution as a function of space for fixed time.

thumbnailFigure 2. Drawdown as a function of space and time.

thumbnailFigure 3. Drawdown as a function of time for a fixed distance.

thumbnailFigure 4. Drawdown as a function of space for fixed time.

7 Conclusion

In this paper, the groundwater flow equation was generalized using the concept of a variational order derivative. The new equation was solved numerically via the Crank-Nicholson technique. We presented in detail the stability and the convergence of this problem. We presented numerical simulations.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

AA wrote the first draft and FJB corrected the final version. Both authors read and approved the final draft.

Acknowledgements

Dedicated to Professor Hari M Srivastava.

The authors would like to thank the referee for some valuable comments and helpful suggestions. This study was supported by the National Research Fund of South Africa.

References

  1. Botha, JF, Buys, J, Verwey, JP, Tredoux, G, Moodie, JW, Hodgkiss, M: Modelling groundwater contamination in the Atlantis aquifer. WRC Report No 175/1/90. Water Research Commission, P.O. Box 824, Pretoria 0001 (2004)

  2. Van Tonder, GJ, Botha, JF, Chiang, WH, Kunstmann, H, Xu, Y: Estimation of sustainable yields of boreholes in fractured rock formations. J. Hydrol.. 241, 70–90 (2001). Publisher Full Text OpenURL

  3. Barker, JA: A generalised radial flow model for hydraulic tests in fractured rock. Water Resour. Res.. 24(10), 1796–1804 (1988). Publisher Full Text OpenURL

  4. Van Der Voort, I: Risk-based decision tools for managing and protecting groundwater resources. PhD dissertation, University of the Free State, P.O. Box 339, Bloemfontein (2001)

  5. Atangana, A: Numerical solution of space-time fractional order derivative of groundwater flow equation. Istanbul, June 20-24. (2012)

  6. Botha, JF, Cloot, AH: A generalized groundwater flow equation using the concept of non-integer order. Water SA. 32(1), 1–7 (2006)

  7. Solomon, TH, Weeks, ER, Swinney, HL: Observation of anomalous diffusion and Lévy flights in a two-dimensional rotating flow. Phys. Rev. Lett.. 71, 3975–3978 (1993). PubMed Abstract | Publisher Full Text OpenURL

  8. Magin, RL, Abdullah, O, Baleanu, D, Zhou, XJ: Anomalous diffusion expressed through fractional order differential operators in the Bloch-Torrey equation. J. Magn. Reson.. 190, 255–270 (2008). PubMed Abstract | Publisher Full Text OpenURL

  9. Sun, HG, Chen, W, Chen, YQ: Variable order fractional differential operators in anomalous diffusion modeling. Phys. A. 388, 4586–4592 (2009). Publisher Full Text OpenURL

  10. Umarov, S, Steinberg, S: Variable order differential equations and diffusion with changing modes. Z. Anal. Anwend.. 28, 431–450 (2009)

  11. Ross, B, Samko, SG: Fractional integration operator of variable order in the Hölder space <a onClick="popup('http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M90','MathML',630,470);return false;" target="_blank" href="http://www.boundaryvalueproblems.com/content/2013/1/53/mathml/M90">View MathML</a>. Int. J. Math. Math. Sci.. 18, 777–788 (1995). Publisher Full Text OpenURL

  12. Chen, YQ, Moore, KL: Discretization schemes for fractional-order differentiators and integrators. IEEE Trans. Circ. Syst. I: Fund. Theoret. Appl.. 49, 363–367 (2002). Publisher Full Text OpenURL

  13. Theis, CV: The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage. Trans. Am. Geophys. Union. 16, 519–524 (1935). Publisher Full Text OpenURL

  14. Zhang, Y: A finite difference method for fractional partial differential equation. Appl. Math. Comput.. 215, 524–529 (2009). Publisher Full Text OpenURL

  15. Tadjeran, C, Meerschaert, MM, Scheffler, HP: A second order accurate numerical approximation for the fractional diffusion equation. J. Comput. Phys.. 213, 205–213 (2006). Publisher Full Text OpenURL

  16. Meerschaert, MM, Tadjeran, C: Finite difference approximations for fractional advection dispersion equations. J. Comput. Appl. Math.. 172, 65–77 (2004). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  17. Deng, WH: Numerical algorithm for the time fractional Fokker-Planck equation. J. Comput. Phys.. 227, 1510–1522 (2007). Publisher Full Text OpenURL

  18. Li, CP, Chen, A, Ye, JJ: Numerical approaches to fractional calculus and fractional ordinary differential equation. J. Comput. Phys.. 230, 3352–3368 (2011). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  19. Chen, CM, Liu, F, Turner, I, Anh, V: A Fourier method for the fractional diffusion equation describing sub-diffusion. J. Comput. Phys.. 227, 886–897 (2007). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  20. Yuste, SB, Acedo, L: An explicit finite difference method and a new Von Neumann-type stability analysis for fractional diffusion equations. SIAM J. Numer. Anal.. 42, 1862–1874 (2005). Publisher Full Text OpenURL

  21. Podlubny, I, Chechkin, A, Skovranek, T, Chen, YQ, Vinagre Jara, BM: Matrix approach to discrete fractional calculus II: partial fractional differential equations. J. Comput. Phys.. 228, 3137–3153 (2009). Publisher Full Text OpenURL

  22. Hanert, E: On the numerical solution of space time fractional diffusion models. Comput. Fluids. 46, 33–39 (2011). Publisher Full Text OpenURL

  23. Lin, R, Liu, F, Anh, V, Turner, I: Stability and convergence of a new explicit finite-difference approximation for the variable-order nonlinear fractional diffusion equation. Appl. Math. Comput.. 212, 435–445 (2009). Publisher Full Text OpenURL

  24. Crank, J, Nicolson, P: A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type. Proc. Camb. Philol. Soc.. 43(1), 50–67 (1947). Publisher Full Text OpenURL

  25. Diethelm, K, Ford, NJ, Freed, AD: Detailed error analysis for a fractional Adams method. Numer. Algorithms. 36, 31–52 (2004)

  26. Li, CP, Tao, CX: On the fractional Adams method. Comput. Math. Appl.. 58, 1573–1588 (2009). Publisher Full Text OpenURL