A generalised groundwater flow equation using the concept of non-integer order derivatives

The classical Darcy law is generalised by regarding the water flow as a function of a non-integer order derivative of the piezometric head. This generalised law and the law of conservation of mass are then used to derive a new equation for groundwater flow. Numerical solutions of this equation for various fractional orders of the derivatives are compared with experimental data and the Barker generalised radial flow model for which a fractal dimension for the flow is assumed.


Introduction
A problem that arises naturally in groundwater investigations is to choose an appropriate geometry for the geological system in which the flow occurs. For example, one can use a model based on percolation theory to simulate the flow in a fractured rock system with a very large fracture density (Berkowitz and Balberg, 1993) or the parallel plate model (De Marsily, 1986) to simulate flow through a single fracture. However, there are many fractured rock aquifers where the flow of groundwater does not fit conventional geometries (Black et al., 1986). This is in particular the case with the Karoo aquifers in South Africa, characterised by the presence of a very few bedding parallel fractures that serve as the main conduits of water in the aquifers (Botha et al., 1998). Attempts to fit a conventional radial flow model to the observed drawdown, see (Van Tonder et al., 2001), always yield a fit that underestimates the observed drawdown at early times and overestimates it at later times, as illustrated in Fig. 1.
The deviation of observations from theoretically expected values is usually an indication that the theory is not implemented correctly, or does not fit the observations. To investigate the first possibility Botha (Botha et al., 1998) 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: (1) where: S 0 is the specific storativity K the hydraulic conductivity tensor of the aquifer Ф(x,t) the piezometric head f(x,t) the strength of any sources or sinks, with x and t the usual spatial and time coordinates V the gradient operator ∂ t 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 (Van Tonder et al., 2001) suggest that the flow is also influenced by the geometry of the bedding parallel fractures, a feature that equation (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 (Barker, 1988) 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 (Van der Voort, 2001), it introduces parameters for which no sound definition exist in the case of non-integer flow dimensions.
As a review of the derivation of Eq. (1) is used as a keystone in the derivation of Eq. (1). This law, proposed by Darcy early in the 19 th 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. 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 idealisation of the reality. Nevertheless, the relative success achieved by Botha, (Botha et al., 1998), 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 reduce to the classical form under the more common conditions. Because K is essentially determined by the permeability of the rocks, and not the flow pattern, the gradient term in Eq.
(2) is the most likely cause for the deviation between the observed and theoretical drawdown observed in the Karoo formations. In this work this possibility is further investigated for a radially symmetric form of Eq. (1), by replacing the classical first order derivative of the piezometric head by a complementary derivative. 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.

Fractional order derivatives
The concept of fractional order derivatives for a function, say f(x), is based on a generalisation of the Abel integral (Courant and John, 1974): ( where: n is a non-zero positive integer Γ(n) the Gamma function (Arfken, 1985) This represents an integral of order n for the continuous function f(x), whenever f and all its derivatives vanish at the origin, x = 0. This result can be extended to the concept of an integral of arbitrary order l, defined as: (4) where: l is a positive real number p an integer such that 0 < r ≤ 1.
Let m now be the least positive integer larger than μ such that μ = m -ρ; 0 < ρ ≤ 1. Equation (4) can then be used to define the derivative of (positive) fractional order, say μ, of a function f(x) as: or equivalently, after performing integration by parts: Note that these results, like Abel's integral, are only valid subject to the condition that f (k) (x)| x=0 = 0 for k = 0,1,2,…,m.
The concept of fractional derivatives is illustrated in Table 1, which lists a few fractional derivatives of the function: The constraint that f(x) and all its derivatives up to the order m must vanish at the origin for D m-ρ f(x) to exist, restricts however the practical application of fractional derivatives considerably. For example, 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. In such situations where the distribution of the piezometric head in the aquifer is a decreasing function of the distance from the borehole the problem may be circumvented by rather using the complementary, or Weyl, fractional order derivative: (7) or, integrating Eq. (7) by parts: where: f(x) satisfies the boundary conditions , (k = 0,1,......,m) It is clear from this definition that when the order of the derivative reduces to an integer, m-1, one has: The way that the fractional and complementary fractional derivatives interpolate the derivatives between integer orders is illustrated in Fig. 1 for the function: As defined above, the fractional derivatives only exist when 0 < ρ ≤ 1. One could in principle also define it for ρ = 0, through the principal Cauchy value of the integral: However, this is not necessary, because the m th fractional derivative with ρ = 0 is identical to the m th fractional derivative in which m is replaced by m + 1 and ρ= 0 by ρ = 1. For example:

A generalised mathematical groundwater flow model
For the sake of clarity the generalisation of the classical model for groundwater flow in the case of density-independent flow in a uniform homogeneous aquifer (Eq. (12)) is considered in this work: where both the specific storativity, S 0 , and hydraulic conductivity, K, are scalar and constant quantities and n=1,2 or 3 is the radial dimension. To be complete the following set of initial and boundary conditions is added: where: Q is the discharge rate of a borehole, with radius r w d the thickness of the aquifer from which the borehole taps water Φ 0 is the initial piezometric head in the aquifer.
In order to include explicitly the possible effect of flow geometry into the mathematical model the radial component of the gradi-ent of piezometric head, ∂ r Φ is replaced by the complementary fractional derivatives of order μ = m -ρ,: This provides a generalised form of the classical equation governing the flow (Eq. (12)): This integro-differential equation does contain the additional parameter μ = m -ρ, 0<ρ≤ 1, which can be viewed as a new physical parameter that characterises the flow through the geological formations. The same transformation generates also a more general form for the boundary condition at the borehole: Relations (14) and (15), together with the initial condition described in Eq. (13), represent a complete set of equations for which a solution exists. The integro-differential character of the relations makes the search for analytical solutions for the problem very difficult however. Therefore this system of equations is solved numerically. A numerical procedure specially developed for this purpose is presented in the next section.

Solution of the generalised flow equation Numerical approximation
The generalised flow (Eq. (14)), subject to the initial and boundary conditions in Eqs. (13) and (15), does not have, to the knowledge of the authors, a known analytical solution. The following two-step numerical approximation is consequently used to study its solution. The first step is to discretise the range in the integral: into a set of line elements ∆r j , j = 1,2,… with equal sizes, while Φ(s, t) is approximated over the element (r j , r j+1 ) with its average Φ (s,t) to obtain: where: Φ l (t) = Φ(r l ,t).
In a second step this approximation and the implicit backward finite difference approximation for the time derivative are used to approximate Eq. (14) as: For example, in the case where n = 1, m = 1 and 0 < ρ ≤ 1, the first order forward difference approximation of the factor V is given by: which reduces in the case where ρ = 1, i.e. when all the coefficients c l are all equal to one, to: the second order forward difference approximation of . On the other hand, if ρ → 0 the c l are all equal to zero and V assumes the form: that may be viewed as a finite difference approximation of over the interval [r i , r i+3 ].

General results
The numerical approximation in the previous section and the parameters summarised in Table 2 are used to compare the solution of Eq. (14) with that of Eq. (12), with results shown in Fig. 3. When pumping commences the perturbation of pressure develops sharply in the vicinity of the borehole and the amplitude of the drawdown is a decreasing function of the order of the derivative considered for the generalised Darcy's law i.e., the maximum value of computed drawdown is obtained when using a zero-order derivative, m= ρ = 1, in the Darcy law.
In that case, the general flow equation reduces to the hyperbolic differential equation: with the shock wave solution: where: H[.] is the Heaviside function.
In the case of one-and two-dimensional flow, n=1,2, the solutions over long times are found to remain unsteady whatever the order of the derivative considered in this work (Fig. 4). However, in the case of tri-dimensional flow (n=3), steady solutions are reached whenever the order of the derivative is less or equal to unity (0< μ ≤ 1) but are unsteady otherwise (Fig. 5). This result may be seen as a generalisation of the well-known result for the classical radial flow equation (i.e. using a classical Darcy's law) where no (non-trivial) steady solution exists, except in the tri-dimensional case.

Case Study: borehole UOat the campus test site
The Campus Test Site of the University of the Free State, Bloemfontein, South Africa, is underlain by a series of mudstones and sandstones from the Adelaide Subgroup of the Beaufort Group of formations in the Karoo Sequence. There are three aquifers present on the site. The top one, a phreatic aquifer, occurs within the upper mudstone layers on the site. This aquifer is separated from the middle and main aquifer, which occurs in a sandstone layer about 10 m thick, by a layer of carbonaceous shale with thickness of 0.5 to 4 m. The bottom aquifer occurs in the mudstone layers more than 100 m thick, which underlies the sandstone unit. A major characteristic of the main aquifer is the presence of a horizontal fracture that coincides approximately with the centre plane of the middle aquifer and which intersects borehole UO5. The fracture itself has an aperture of approximately 10 mm but the adjacent 200 mm of the sandstone is also highly permeable. In a recent article Van Tonder (Van Tonder et al., 2001) reports data for the drawdown obtained with a pressure transducer during a constant rate test in which UO5 is pumped at a constant rate of 1.25 ℓ·s -1 . An analysis of the data using a Theis type of curve yielded a value of 19 m 2· d -1 for the transmissivity, T, and 8.6 for the storativity, S. This best fitted Theis curve is shown in Fig. 1. It is clear from this graph that the solution of the Theis equation is reasonably matching the amplitude of the measured drawdown locally in time, for several hours after pumping started. However, it is also obvious that this fitting is realised with a function that overestimates the drawdown during the early times and underestimates it at large times. This suggests that the differential operator used by the Theis model does not accommodate a time-dependent solution that can account properly for the dynamical evolution of the piezometric value. One way to circumvent this problem is to enlarge the spectrum of possible solutions supported by the differential operator. This possibility was investigated by re-computing the drawdown using the same storativity and transmissivity value as obtained with the Theis model but assuming that Darcy law is represented by a complementary fractional order derivative, the order of the derivative, as well as the flow thickness, being considered free parameters. The result of this exercise using a order the derivative of μ = 1.15 and a flow thickness d = 0.04 m is shown in Fig. 1.

Comparison with Barker's Model
Based on the assumption that groundwater in an aquifer flows through equipotential surfaces that are projections of l-dimen- sional spheres onto two-dimensional space, Barker, (Barker, 1988), obtains the following analytical solution of Eq.(12) for an infinite aquifer with a line source: where: Γ(.,.) is the incomplete gamma function, λ the dimension of the flow, which equals the spatial dimension n when being an integer .
The other quantities all have the same meanings as before.
Although Barker's model differs conceptually from the one based on the generalised Darcy law introduced above, a closer inspection of the two models shows that they have much in common. For example, both models are based on Eq. (12) and the boundary conditions in Eq. (13). The real differences lie in the interpretation of the parameters λ and n and the flow areas, which assumes respectively the forms: Beside the fact that the two models yield identical results whenever λ = n and m =2; ρ = 1, they also possess quite similar solutions in the cases where non-integer values of the parameters are considered. As for example, in the case of bi-dimensional (or a tri-dimensional) flows, there exists an inverse relationship between the order of derivatives μ used in the generalised Darcy model and the flow dimension λ considered in Barker's (see Fig. 6). An interesting feature is that this relationship is not restricted to order of derivatives larger than one but is also observed for values ranging in 0< μ < 1. In that case it must be assumed that a dimension higher than 3 in Barker's model is feasible.

Discussion
The expression of the Darcy velocity as a function of a complementary derivative of fractional order may sound to be merely some kind of highly technical mathematical concept and before drawing conclusions on this new model, one would like to attract the attention of the reader on some basic physical features implied by the model. Broadly seen the generalised version of the Darcy law is affecting conceptually the physical content of Darcy's work in mainly two ways. Firstly, the new generalised constitutive equation is not, unlike the usual Darcy relation, satisfying the popular principle of local (spatial) action that postulates that the evolution of a system at a given spatial position and at a given time is only affected by the behaviour of the different variables in the direct neighbourhood around that point. Within the frame of Darcy's law this principle implies that the Darcy velocities may only be expressed as relations containing the value, at that point, of the piezometric head and/or its derivatives. The use of a constitutive relation based on a spatial integration of the piezometric head over the half infinite line is breaking away from this postulate and is allowing for the consideration of more general situations where the groundwater flow at a given point of the aquifer is governed not only by the properties of the piezometric field at that specific position but also depends on the global spatial distribution of that field in soil matrix. It is however expected that the information relative to , ( 4 7 the behaviour of the fluid particle presently at the sampling point for different times into the past, t < γ. In that case, the collection of information at an increasingly large distance is equivalent to taking a longer time period of the fluid particle history into account. Of course, unless the global behaviour of the groundwater flow is steady, those images are blurred pictures of the past of the fluid particle and it may be expected that images of the far past are more distorted than those relative to a more contemporary history. For that matter the weighting factor contained in the integral model (Eq. (15)) is also acting as a time filter according to which only the recent past (i.e. the less blurred) is taken into account. In the case of a steady motion, the trajectories followed by the fluid particles with respect to time identify with the curves made of the fluid particles that, at a given time, will go or went through a given point in the space and therefore the behaviour of the flow along the half infinite line right to the sampling point is an exact representation of the history of the fluid particle at the sampling point. A direct consequence of the introduction of the notion of history into the generalised model is that the global equation describing the groundwater flow is loosing its parabolic character and the speed at which the information propagates throughout the porous medium is finite, as expected. A second main difference between the classical Darcy law and the generalised model is that the intrinsic geometry of a specific aquifer is taken into account. Generally speaking, the geometry of the medium is accounted for in the model in two different ways. On the one hand, the effect of the geometry of the flow is implicitly included in the model through the integral character of the equation. As explained above the integral model is referring to the "space-like" history of the variables that, in its turn, is accounting for the geometry of the flow as experienced by the fluid particle up to date. Therefore the geometry of the aquifer lying ahead of a given particle has no influence on the behaviour of that particle. From the other hand, the fractional order derivative is accounting explicitly for the geometry of the aquifer by introducing an additional coefficient, in the definition of the Darcy velocities. This new quantity may be combined with the classical hydraulic conductivity, K, to define a more comprehensive generalised hydraulic parameter Kf: It is worth noticing that this geometrical factor is independent of the physical dimension assumed for the flow description (i.e. linear, radial or spherical flow) and that it is always present in the model, except in the special case μ = 1.

Conclusions
The classical Darcy Law has been generalised using the concept of (complementary) fractional order derivatives. This leads to the formulation of a new (generalised) form of the groundwater flow equation. For constant rate test, numerical simulations performed show that, likewise in the case of the classical Darcy Law, a non-trivial steady state solution is only possible when the dimension of the flow is greater than two and the order of the derivative is taken lesser or equal to one. A comparison of these results 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.