Abstract
The classical Darcy law is generalized by regarding the water flow as a function of a noninteger order derivative of the piezometric head. This generalized law and the law of conservation of mass are then used to derive a new equation for groundwater flow. Two methods including Frobenius and Adomian decomposition method are used to obtain an asymptotic analytical solution to the generalized groundwater flow equation. The solution obtained via Frobenius method is valid in the vicinity of the borehole. This solution is in perfect agreement with the data observed from the pumping test performed by the institute for groundwater study on one of their boreholes settled on the test site of the University of the Free State. The test consisted of the pumping of the borehole at the constant discharge rate and monitoring the piezometric head for 350 minutes. Numerical solutions obtained via Adomian method are compared with the Barker generalized radial flow model for which a fractal dimension for the flow is assumed. Proposition for uncertainties in groundwater studies was given.
1. Introduction
The real problem encounter in groundwater studies up to now is the real shape of the geological formation in which water flows in the aquifer under investigation. However, there are many fractured rock aquifers where the flow of groundwater does not fit conventional geometries [1], for example, in South Africa, the Karoo aquifers, characterized by the presence of a very few bedding parallel fractures that serve as the main conduits of water in the aquifers [2]. With a challenge to fit the solution of the groundwater flow equation with experimental data from field observation in particular, the observed drawdown see [3], for all time yields a fit that undervalues the observed drawdown at early times and overvalues it at later times. The variation of observations from theoretically predictable values is usually an indication of uncertainties in the predictable. To investigate the first possibility Botha et al. [2] developed a three-dimensional model for the Karoo aquifer on the campus of the University of the Free State. This model is based on the conventional, saturated groundwater flow equation for density-independent flow: where is the specific storativity, the hydraulic conductivity tensor of the aquifer, the piezometric head, ) the strength of any sources or sinks, with and the usual spatial and time coordinates; the gradient operator, and the time derivative.
This model showed that the dominant flow field in these aquifers is vertical and linear and not horizontal and radial as commonly assumed. However, more recent investigations [3] suggest that the flow is also influenced by the geometry of the bedding parallel fractures, a feature that (1) cannot account for. It is therefore possible that the equation may not be applicable to flow in these fractured aquifers.
In an attempt to circumvent this problem, Barker [4] introduced a 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 noninteger flow dimensions.
As a review of the derivation of (1) will show, [5] Darcy law is used as a keystone in the derivation of (1). This law, proposed by Darcy early in the 19th century, is relying on experimental results obtained from the flow of water through a one-dimensional sand column, the geometry of which differs completely from that of a fracture [6]. There is therefore a possibility that the Darcy law may not be valid for flow in fractured rock formations but is only a very crude idealization of the reality [6]. Nevertheless, the relative success achieved by Botha et al. [2], to describe many of the properties of Karoo aquifers, suggests also that the basic principle underlying this law may be correct: the observed drawdown is to be related to either a variation in the hydraulic conductivity of the aquifer or a change in the piezometric head. Any new form of the law should therefore be reduced to the classical form under the more common conditions. Because is essentially determined by the permeability of the rocks, and not the flow pattern, the gradient term in (2) is the most likely cause for the deviation between the observed and theoretical drawdown observed in the Karoo formations [6]. In the same direction, Cloot and Botha introduced the concept of non-integer fractional derivative to investigate a radially symmetric form of (1), by replacing the classical first order derivative of the piezometric head by a complementary fractional derivative [6]. However the generalized model for groundwater flow equation was solved numerically. In this work, more general form of groundwater flow equation will be introduced; the Frobenius and Adomian decomposition will be used to give an asymptotic solution of the generalized model for groundwater flow equation. Because the concepts of fractional (or non-integer) order derivatives and complementary fractional order derivatives may not be widely known, both concepts are first briefly discussed below.
2. Fractional Order Derivatives
On one hand, the concept of fractional calculus is popularly believed to have stemmed from a question raised in the year 1695 by Marquis de L Hospital (1661–1704) to Gottfried Wilhelm Leibniz (1646–1716), which sought the meaning of Leibniz’s currently popular notation for derivative of order when (what if ). In his reply, dated September 30, 1695, Leibniz wrote to L’Hospital as follows: “This is an apparent paradox from which, one day, useful consequences will drawn.” On the other hand, the concept of fractional order derivatives for a function, say , is based on a generalization of the Abel integral [7]: where is a nonzero positive integer [8–12] and the Gamma function [13]. This represents an integral of order for the continuous function , whenever and all its derivatives vanish at the origin, . This result can be extended to the concept of an integral of arbitrary order , defined as: where is a positive real number and an integer such that .
Let now be the least positive integer larger than such that ; . Equation (4) can then be used to define the derivative of (positive) fractional order, say , of a function as Note that these results, like Abel’s integral, are only valid subject to the condition that for [8–11].
2.1. Properties
Properties of the operator can be found in [14, 15], we mention only the following: for , , and :
3. A Generalized Mathematical Groundwater Flow Model
For the sake of clarity the generalization of the classical model for groundwater flow in the case of density-independent flow in the uniform homogeneous aquifer is considered in this paper [4]. Consider the following groundwater flow equation where both the specific storability, , and hydraulic conductivity, , are scalar and constant quantities and or is the radial dimension. To be complete, the following set of initial and boundary conditions is added: Here means that before a pumping test begins the level of water or the initial hydraulic head in the aquifer is linear function of space with a positive gradient to be found, means that during the pumping test the level of water is not affected for a very long distance from the borehole, and means that the rate of pumping is proportional to the hydraulic conductivity.
Hence is the discharge rate of the borehole, with radius and the thickness of the aquifer from which the borehole taps. In order to include explicitly the possible effect of the geometry into mathematical model the radial component of the gradient of the piezometric head, is replaced by the Weyl-fractional derivatives of order ; the fractional derivative in this paper is Caputo derivative and is defined as follows: This provides a generalized form of the classical equation governing the flow equation (1): The integrodifferential equation does contain the additional parameter , , which can be viewed as a new physical parameter that characterizes the flow through the geological formations [6]. The same transformation generates also a more general form for the boundary condition at the borehole [6]: The relations (10) and (11), together with the initial condition described in (10), represent a complete set of equations for which a solution exists. The integro-differential character of the relations makes the search for analytical solution very difficult however. Nevertheless, in this paper we make use of Frobenius and Adomian decomposition methods to give an asymptotic solution.
4. Solution of the Generalized Flow Equation
4.1. Frobenius Method
In this work to perform the Frobenius method, we consider the groundwater flow governed by the following fractional Caputo-Weyl derivation partial differential of order , where is real number, . Also, we consider the dimension of the flow to be 2. Therefore, for (10) can then be transformed into the following equation: Applying the Laplace operator on both sides of the above equation, we have the following ordinary differential equation where is the variable of Laplace. In this matter we choose or we choose , meaning that the level of water is the same everywhere in the aquifer if the water is not taken out from the aquifer. If we let then (13) becomes We put and . To meet the condition under which Frobenius method can be used, we have to prove that and are analytical around that means we have to prove that and can be written as series. It is very obvious to see that and can be expressed as follows: We start here with the coefficients of : implying that .
Putting we have that and equating the coefficients of same power we obtain the following set of equations: Therefore, the coefficients can be given with the general following recursive formula: And obviously the coefficients are given below as From the above expression we can see that and are analytical around , which follows from Frobenius method that the solution of (15) can be in the form This solution is not convergent for a large ; that is, the solution will diverge if we observe the drawdown at a position very far from the borehole from which the water is taken out. Therefore, we restrict our solution in the vicinity of the borehole; more precisely we investigate the solution when . That means we pump water from the borehole and we observe the drawdown in the vicinity of the borehole.
Thus substituting (21), , and into (15) and equating the coefficients of the same power, we obtain the following recursive formula for which the coefficients , coefficients of our series: Here we have the following: However, making use of the boundary condition in (11), we have the following: Then the general solution of the fractional Caputo-Weyl derivation of groundwater flow of order in Laplace space is given below as: To observe the behavior of the solution at the borehole, which corresponds to , the series solution is reduced here to the coefficient with order zero which is obtained when and it’s given below as In the following section the analytical asymptotic solution obtained in Laplace space via Frobenius method will be compared with the experimental data.
4.2. Numerical Results
In order to examine the validation of this solution, the above asymptotic solution is compared with the experimental data from the pumping test performed by the Institute for Groundwater Studies on one of their borehole settled on the campus test site of the University of the Free State. The test consisted of the pumping of the borehole at the constant discharge rate and monitoring the piezometric head for 350 minutes. The first step in the section is to discretize the range of Laplace transform since the exact Laplace transform cannot be obtained in practice. This is done as follows: For we approximate where for ; then the results obtained in the fields and Laplace transform become Using this numerical scheme, the physical data was transformed into Laplace space. A comparison between these values and asymptotic computed data can only be provided in the in real space not in Laplace space. Since it is not worth concluding the validity of this solution in Laplace space, the inverse Laplace transform is applied in (26). To test the validity of this solution in real space and applying the inverse Laplace transform on (26), (29) is obtained
The above solution is compared graphically to the experimental data from the pumping test performed by the institute for groundwater study on one of their borehole settled on the test site of the University of the Free State. The small difference observed in the above graph (Figure 1) is due to uncertainties in measurement and this will be discussed in Section 6 of this work. The aquifer parameters used in this models are recorded in Table 1, the observed data from field observation will be attached to this paper. Although this solution is in agreement with the experimental data, there will be the need to investigate the case where the observation can be done for a long distance. In the next section another approach will be introduced to solve the space-time-fractional derivative of groundwater flow equation, and this method is the Adomian decomposition method.
5. Adomian Decomposition Method
5.1. Example 1
This section is concerned with the groundwater equation with time- and space-fractional derivatives of the form Subject to the initial and boundary conditions described in (8), the level of water is assumed to be the same throughout the aquifer before the pumping so that the gradient described in (8) is zero. Furthermore, it is assumed that a fractional change in drawdown is constant for meaning .
The method used here is based on applying the operator on both sides of (30) to obtain The Adomian decomposition method [16, 17] assumes a series solution for (31) to be where the components are determined recursively. Substituting (32) into both sides of (30) gives Following the decomposition method, the recursive relations are introduced as It is worth noting that if the component is defined, then the remaining components can be completely determined such that each term is determined by using the previous terms, and the series solutions are thus entirely determined. Finally, the solution is approximated by the truncated series However, the inclusion of boundary conditions in fractional differential equations introduces additional difficulties. The Adomian decomposition method can handle these difficulties by using the space-fractional operator and the initial conditions only. The method provides the solution in the form of a rapidly convergent series that may lead to the exact solution in the case of integer derivatives and to an efficient numerical solution with high accuracy for 0 . The convergence of the decomposition series has been investigated in [18, 19].
Following the recursive formula equation (35) and using the fact that , the equations below are obtained: The component was also determined and will be used, but for brevity it is not listed. In this matter five components of the decomposition series (30) were obtained for which was evaluated to have the following expansion: Applying the boundary condition yields where
5.2. Example 2
Consider the groundwater flow equation with time- and space-fractional derivatives of the form subject to the initial condition described in (11). Furthermore we suppose that the gradient , meaning that the water level is everywhere in the aquifer at and the fractional change in drawdown is a constant and boundary condition: Following the discussion earlier, we have the below recursive formula: It follows from the recursive formula that The component was also determined and will be used, but for brevity it is not listed. In this matter five components of the decomposition series (41) were obtained of which was evaluated to have the following expansion: Applying the boundary condition yields where
5.3. Example 3
Consider time fractional derivative of the groundwater flow equation with time-fractional derivatives of the form Subject to the initial and boundary conditions described in (8) it is assumed that . Following the discussion presented earlier, we obtained the below recursive formula: The component was also determined and will be used, but for brevity not listed. In this matter five components of the decomposition series (31) were obtained of which was evaluated to have the following expansion Applying the boundary conditions yield to where The normalized solutions with of (50), (51) and (52) are illustrated graphically in Figures 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, and 13 for different values of various parameters; these solutions are compared to the solution proposed by Barker [4]. These graphs show the behaviour of the drawdown during the pumping test, first as a function of space and time from the borehole to the point of observation, secondly as function of time for a fixed distance from the borehole, and finally as a function of space for a fixed time. Based on the assumption that groundwater in an aquifer flows through an equipotential surface that are projections of -dimensional spheres onto two-dimensional space, Barker [4] obtained the following analytical solution for an infinite aquifer with a line source: where is the incomplete gamma function and the dimension of the flow which equals the special dimension, being an integer, and the other quantities all have the same meaning as before. 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. The main raison of comparison of these results obtained via the Adomian decomposition methods with those of Barker’s fractal radial flow model is to establish a possible relationship between the fractional order of the derivative and the parameter fractal introduced earlier by Barker.
Table 2 shows the theoretical values of the discharge rate and aquifer parameters used in the numerical simulations.
6. Discussion and Propositions
Although the analytical solution obtained via Frobenius method fit the experimental data or have described successfully the events taking place in the vicinity of the borehole, on one hand, and the analytical solutions obtained via Adomian decomposition method was successfully compared to the solution proposed by Barker, on the other hand, the problem of choosing an appropriate geometry for the geological system in which the flow occurs still remains a challenge in groundwater studies. We personally believe that we do not propose solution to a problem because it going to be useful for the generation in which we are, but we propose solution because we hope that they will one day be useful for mankind; therefore the following proposition may not be useful for this generation because we may not have the adequate technology to perform the steps involved but it will be in the future. We think that describing the groundwater flow with one equation for the whole aquifer is unrealistic, because from one point of the aquifer to another properties including geology and geometry change, therefore the flow. Assumption such as homogeneity, isotropy, uniform thickness over the area under investigation and so on render the study of groundwater uncertain. In order to study the geology or the geometry of an aquifer, we have to divide the aquifer in small portions, from north to south and from east to west. The geology and the geometry of each portion can then be studied. If the results of the study reveal that the portion under investigation is for instance a hyperboloid, since there is no exact solution to groundwater flow model by hyperboloid flow and there are solutions to circular flow, then a suitable transformation can be done including transformation of hyperboloid coordinates to Cartesian coordinates, then to Cartesian coordinates to cylindrical coordinate such that this new coordinate can now be put in the equation describing groundwater flow, and the solution to the new equation can then be investigated. Henceforth knowing the real geology and geometry of this portion, the real paths flow will be known. Then having a good knowledge of each small portion including its geometry and geology, the real geometry and geology of the aquifer can be not exactly but more accurately known.
For groundwater remediation it will be possible to know where the maximum, minimum, and average chemical concentrations are found in the aquifer. We have no proof for this but we believe that this proposition will be useful in reducing uncertainties in groundwater study. It is believed that the field test gives the characteristic of an aquifer, but we believe that the field test gives both uncertainties and characteristics of the aquifer; therefore quantify uncertainties in this measurement lead us to the real picture of aquifer characteristics, henceforth we propose that the studies in groundwater should focus on both uncertainties and fields observations, because what is known is bounded by what is not known; knowing what is not known give a real picture to what was known, and it follows that the knowledge of uncertainties in groundwater study will give a clear picture of what we already know in groundwater.
7. Conclusion
The classical Darcy Law has been generalized using the concept of complementary fractional order derivatives of Weyl fractional derivative. This leads to the formulation of a new generalized form of the groundwater flow equation [6]. The applications of Adomian decomposition and Frobenius methods were extended to obtain explicit and numerical solutions of the space-time fractional groundwater flow. The two methods were very clearly efficient and powerful techniques in finding the solutions of the proposed equations. The solution obtained via Frobenius takes into account the events taking place in the vicinity of the borehole during the pumping test whereas the solution obtained via Adomian decomposition methods takes into account the events that take place far from the point where water is pumped out, that is, in the borehole. The solution obtained via Frobenius method was in perfect agreement with the observed data obtained from the pumping test performed by the Institute for Groundwater Studies on one of their borehole settled on the campus test site of the University of the Free State. The Adomian decomposition method requires less computational work than existing approaches while supplying quantitatively reliable results. The obtained results demonstrate the reliability of the algorithms and their wider applicability to fractional evolution equations. A comparison of these results obtained via the Adomian decomposition methods with those of Barker’s fractal radial flow model suggests that there exists a relation between the fractional order of the derivative and the non-integral dimension of the flow.
Authors’ Contribution
Abdon Atangana wrote the first draft and P. D. Vermeulen read and revised it; the revised version was read and corrected by both authors.
Conflict of Interests
The authors declare that there is no conflict of interests for this paper.