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
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  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:
Figure 1. Groundwater flow in the borehole during pumping test.
Here is the specific storativity, is the hydraulic conductivity tensor of the aquifer, Φ is the piezometric head, is the strength of any sources or sinks. As a review of the derivation of Eq. (1.1), we will show the Darcy law
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 has the dimensions of velocity. Recent investigations  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  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 , 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 .
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 . 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 , denote a continuous and necessary differentiable, let be a continuous function in . Then its variational order differential in is defined as
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)  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 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
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 . Yuste et al. introduced weighted average finite difference methods . Podlubny et al. proposed the matrix approach for fractional diffusion equations , and Hanert proposed a flexible numerical scheme for the discretization of the space-time fractional diffusion equation . Recently, Zhuang et al. considered the numerical schemes for VO space fractional advection-dispersion equation , Lin et al. investigated the explicit scheme for VO nonlinear space fractional diffusion equation . 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 , , , , , , h is the step and τ is the time size, M and N are grid points.
3.1 Crank-Nicholson scheme 
We introduce the Crank-Nicholson scheme as follows. Firstly, the discretization of first- and second-order space derivative is stated as follows:
The Crank-Nicholson scheme for the VO time fractional groundwater flow model can be stated as follows:
Now, replacing equations (3.1), (3.2), (3.3) and (3.4) in (2.3), we obtain the following:
For simplicity, let us put:
Equation (3.5) becomes
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 . Here 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 Chen et al. that
Observe that for all , , in addition, according to the problem in point, the storativity S, and transmissivity T are positive constants. Then the following properties of the coefficients , , , and can be established:
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 . 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:
If we assume that in equation (4.1) can be put in the delta-exponential form as follows:
where ϕ is a real spatial wave number, new replacing the above equation (4.5b) into (4.5a), we obtain
Equation (4.6) can be written in the following form:
Our next concern here is to show that for all , the solution of equation (4.7) satisfies the following condition:
To achieve this, we make use of the recurrence technique on the natural number k.
For and remembering that , are positive for all , we obtain
Assuming that for the property is verified, we get
Making use of the triangular inequality we obtain
Using the recurrence hypothesis, we have
which this completes the proof.
5 Convergence analysis of the Crank-Nicholson scheme
If we assume that ( , ) is the exact solution of our problem at the point , by letting and , substituting this in equation (4.7), we obtain
From equations (3.1) and (3.4), we have
From the above, we have that
where , , 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. and further work by Li and Tao .
Lemma 1 is true for ( ), where , Kis a constant. In addition,
This can be achieved via the recurrence technique on the natural numberk.
When , we have the following:
Now suppose that , . Then
which completes the proof.
Theorem 1The Crank-Nicholson scheme is convergent, and there exists a positive constantVsuch that
An interested reader can find the solvability of the Crank-Nicholson scheme in the work done by . 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 , , , and .
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.
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.
The authors declare that they have no competing interests.
AA wrote the first draft and FJB corrected the final version. Both authors read and approved the final draft.
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.
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)
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
Barker, JA: A generalised radial flow model for hydraulic tests in fractured rock. Water Resour. Res.. 24(10), 1796–1804 (1988). Publisher Full Text
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
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
Sun, HG, Chen, W, Chen, YQ: Variable order fractional differential operators in anomalous diffusion modeling. Phys. A. 388, 4586–4592 (2009). Publisher Full Text
Ross, B, Samko, SG: Fractional integration operator of variable order in the Hölder space . Int. J. Math. Math. Sci.. 18, 777–788 (1995). Publisher Full Text
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
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
Zhang, Y: A finite difference method for fractional partial differential equation. Appl. Math. Comput.. 215, 524–529 (2009). Publisher Full Text
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
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
Deng, WH: Numerical algorithm for the time fractional Fokker-Planck equation. J. Comput. Phys.. 227, 1510–1522 (2007). Publisher Full Text
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
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
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
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
Hanert, E: On the numerical solution of space time fractional diffusion models. Comput. Fluids. 46, 33–39 (2011). Publisher Full Text
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
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
Li, CP, Tao, CX: On the fractional Adams method. Comput. Math. Appl.. 58, 1573–1588 (2009). Publisher Full Text