# The Use of Fractional Order Derivative to  Predict the Groundwater Flow

Authors: Abdon Atangana; Necdet Bildik
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

The aim of this work was to convert the Thiem and the Theis groundwater flow equation to the time-fractional groundwater flow model. We first derived the analytical solution of the Theim time-fractional groundwater flow equation in terms of the generalized Wright function. We presented some properties of the Laplace-Carson transform. We derived the analytical solution of the Theis-time-fractional groundwater flow equation (TFGFE) via the Laplace-Carson transform method. We introduced the generalized exponential integral, as solution of the TFGFE. This solution is in perfect agreement with the data observed from the pumping test performed by the Institute for Groundwater Study on one of its borehole 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 rateQ and monitoring the piezometric head for 350 minutes.


---

## Body

## 1. Introduction

Groundwater problem is perhaps one of the most difficult real-world problems to be modelled into mathematical formulation.  To model this problem accurately, one must know precisely the behavior of the medium through which the water is moving. However, this medium through which the flow occurs can change from one point to another, also from one period to another. For example, the hydraulic conductivity of an aquifer can differ from one direction to another. Several scholars have intensively tried to propose a better model that can be used to predict the movement of water through the aquifer. However, their results still present some lacks.  Recently, It was revealed that real problems modelled via fractional order derivative present better results when matching their mathematical representation with experimental data.  To test this, Botha and Cloot [1] presented some good results by generalizing the groundwater flow equation to the concept of fractional order derivatives. In the same line of ideas, Atangana [2] examined an approximate solution of the generalized groundwater flow equation via the Frobenius method. The results obtained from his investigation showed better prediction. Recently, Atangana and Botha further extended the fractional groundwater equation to the concept of the fractional-variation order groundwater [3]. They presented the stability and the convergence of the numerical scheme via Crank-Nicolson method. Up to now, there is no approximate or exact analytical mathematical expression that can be used to describe the solution of the fractional groundwater flow equation. Therefore, one of the purposes of this work is to present some analytical mathematical expression than can be used as approximate solution of the time-fractional groundwater flow equation.An aquifer test (or a pumping test) is conducted to evaluate an aquifer by “stimulating” the aquifer through constant pumping and observing the aquifer's “response” (drawdown) in observation wells. Aquifer testing is a common tool that hydrogeologists use to characterize a system of aquifers, aquitards, and flow system boundaries. Aquifer tests are typically interpreted by using an analytical model of aquifer flow (the most fundamental being the Theis solution) to match the data observed in the real-world then assuming that the parameters from the idealized model apply to the real-world aquifer. In more complex cases, a numerical model may be used to analyze the results of an aquifer test, but adding complexity does not ensure better results. For the most part frequently, an aquifer assessment is carried out by propelling water out from one borehole at a fixed speed and for as a minimum of  24 hours at the same time as cautiously evaluating the water levels in the observed borehole. When water is pumped from the pumping well, the pressure in the aquifer that feeds that well declines. This decline in pressure will show up as drawdown (change in hydraulic head) in an observation well. Drawdown decreases with radial distance from the pumping well and drawdown increases with the length of time that the pumping continues. The aquifer characteristics which are evaluated by most aquifer tests are [4, 5] as follows.(i)
The hydraulic conductivity is the rate of flow of water through a unit crosses sectional area of an aquifer at a unit hydraulic gradient. In English units the rate of flow is in gallons per day per square foot of cross sectional area.(ii)
Specific storage or storativity being a measure of the amount of water a confined aquifer will give up for a certain change in head.(iii)
The transmissivity is the rate at which water is transmitted through a unit thickness of an aquifer under a unit hydraulic gradient. It is equal to the hydraulic conductivity times the thickness of an aquifer.The rest of this paper has been presented as follows. In Section 1, we presented the background of the fractional order derivative. We derived the analytical solution of the Theim fractional groundwater flow equation in Section 3.  In Section 4, we presented the derivation of the Theis fractional groundwater flow via the Laplace-Carson transform method. We presented an alternative analytical solution of the Theis fractional groundwater flow equation in terms of the generalized exponential integral in Section 5. The numerical comparisons with experimental data are presented in Section 6 and the conclusion is in Section 7. We will start with the background of the fractional derivative.


## 2. Background of the Fractional Order Derivative

There exists the vast literature on different definitions of fractional derivatives. The most popular ones are the Riemann–Liouville and the Caputo derivatives [6–11]. Caputo’s definition has the form of
(1)Dxα(f(x))=1Γ(n-α)dndxn×∫0x(x-t)n-α-1f(t)dt,n-1≤α≤n.
For the case of the Caputo fractional order derivative, we have the following definition:
(2)Dxα(f(x))=1Γ(n-α)×∫0x(x-t)n-α-1dndxnf(t)dt,n-1≤α≤n.
Each of the pervious fractional order derivatives presents some advantages and disadvantages [6–10]. The Riemann-Liouville derivative of a constant is not zero while Caputo derivative of a constant is zero but demands higher conditions of regularity for differentiability [6–10]: to compute the fractional derivative of a function in the Caputo sense, we must first calculate its derivative [11]. Caputo derivatives are defined only for differentiable functions while functions that have no first-order derivative might have fractional derivatives of all orders less than one in the Riemann-Liouville sense [12, 13]. Guy Jumarie has recently modified the Riemann-Liouville derivative (see [14])
(3)Dxα(f(x))=1Γ(n-α)dndxn×∫0x(x-t)n-α-1{f(t)-f(0)}dt.
The Caputo fractional derivative will be considered in this work due to the applicability of the Caputo derivative in real-world problems [15].


## 3. Thiem’ Groundwater Flow Equation

The Theim fractional groundwater flow equation is an ordinary differential equation given in the following. The equation describes the change in level of water as function of distance during the pumping test [5]:
(4)DrrαΦ(r)+1rΦ(r)=0,1<α≤2.
Subject to the initial condition, Q=2πTDr(Φ(rb)).We will make use of the Laplace transform to derive analytical solution of (4). Thus, multiplying on both sides of (4) by r and secondly applying the Laplace transform, we obtain the following expression:
(5)d[L(Φ)(s)]ds+(αs+1sα)(L(Φ)(s))=∑m=2ldm(m-1)sm-2-α,
where  dm=D0+α-mΦ(0+)(m=2,…,l). Now, one can derive the solution of the ordinary order differential equation with respect to the Laplace transform of Θ(s)=L(Φ(r)):
(6)Θ(s)=s-αexp[-s1-α1-α]×[a1+∑m=2ldm(m-1)∫sm-2exp[-s1-α1-α]ds],
with a1 an arbitrary real constant that will be obtained via the initial condition. We next expand the exponential function in the integrand in a series, and using term-by-term integration, we arrive at the following expression:
(7)Θ(s)=cΘ1(s)+∑m=2ldm(m-1)Θm*(s)
with, of course,
(8)Θ1(s)=s-αexp[-s1-α1-α],Θm*(s)=s-αexp[s1-αα-1]∑j=0∞(11-α)js(1-α)j+m-1[(1-α)j+m-1]j!.
Now, applying the inverse Laplace transform on Θ1(s) and using the fact that
(9)s-[α+(α-1)j]=L[rα+(α-1)j-1Γ(α+(α-1)j)]
we obtain
(10)Φ1(r)=rα-1oΨ1[(α,α-1)∣xα-1α-1]
with oΨ1[·] the generalized Wright function [16] for p=1 and q=2. We next expand the exponential function exp[-s1-α/(1-α)] in power series; multiplying the resulting two series; in addition of this if we consider the number bk(α,m) defined for α>0, m=2,…l, (α≠(p+m-1)/p,p∉ℕ) and k∈ℕ0,
(11)bk(α,m)=∑p,j=0,…k,p+j=kl(-1)qp!j!(1-α)q+m-1.
The previous family of number possesses satisfies the following recursive formula:
(12)bk(α,m)bk+1(α,m)=α-mα-1+k,
which produces the explicit expression for bk(α,m) in the form of
(13)bk(α,m)=Γ[(α-m)/(α-1)](m-1)Γ[((α-m)/(α-1))+k],k∈ℕ0.
Now, having the previous expression on hand, we can derive that
(14)Θm*(s)=sm-α-1(∑j=0∞(11-α)js(1-α)pp!)×(∑p=0∞(11-α)p(-1)p[(1-α)p+m-1]s(1-α)jp!)=∑k=0∞bk(α,m)(11-α)ks(1-α)k+m-α-1(m=2,…,l).
However, remembering (9) with β=(α-1)k+α+1-m, we can further derive the following expression for Φm*(r) as
(15)Φm*(r)=∑k=0∞bk(α,m)(11-α)k×Γ(k+1)Γ[α+1-m+(α-1)k]x(α-1)k+α-mk!
or in the simplified version we have
(16)Φm*(r)=Γ[(α-m)/(α-1)](m-1)Φm(r),
where
(17)Φm(r)=rα-mΨ21×[(1,1)(α+1-m,α-1),(α-mα-1,1)∣rα-1α-1].
It follows that the solution of the fractional-Thiem groundwater flow equation is in the form of
(18)Φ(r)=a1rα-1oΨ1[(α,α-1)∣xα-1α-1]+a2∑m=2lbm(m-1)Γ[(α-m)/(α-1)](m-1)rα-m×Ψ21[(1,1)(α+1-m,α-1),(α-mα-1,1)∣rα-1α-1].
In our case l=2.


## 4. Time-Fractional Theis Groundwater Flow Equation

The easiest sweeping statement of subsurface water flow equation, which while we are on the subject is in addition in harmony in the midst of the real physics of the observed fact, is to presume that water level is not in a balanced state but momentary state. Theis (1935) [17] 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 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:
(19)SDtΦ(r,t)=TDrrΦ(r,t)+1rDrΦ(r,t).
The aforementioned equation is classified under parabolic equation. To include explicitly the variability of the medium through which the flow takes place, the standard version of the partial derivative respect to time is replaced here with time-fractional order derivative to obtain
(20)SDtαΦ(r,t)=TDrrΦ(r,t)+1rDrΦ(r,t),0<α≤1
with initial condition Φ(r,0)=0 and boundary condition limr→∞Φ(r,t)=0, Q=2πT∂rΦ(rb,t), here T is the transmissivity of the aquifer, rb is the ratio of the borehole, and Q is the discharge rate, or the rate at which the water is being taken out of the aquifer.Some few integral transform operators have been intensively used to solve some kind of ordinary and partial differential equations.  See, for instance, the Fourier transform, the Laplace transform, the Mellin transform [18], and the Sumudu transform [19–22]. Beside these integral operators, there exists a similar operators called the Laplace-Carson transform [23]; this operator has been neglected. However, this operator has some properties that can be used to solve a kind of ODE, PDE, FODE, and FPDE. The aim of this section is therefore devoted to the discussion underpinning the definition, properties of the Laplace-Carson transform, and its application to the fractional groundwater flow equation. We shall start with the definition and properties.Definition 1.
Letf(x) be a continuous function over an open interval (0,∞) such that its Laplace transform is n time differentiable; then the Laplace-Carson transform f is defined as follows:
(21)Lc(s)=Lc[f(x)](s)=∫0∞xe-xsf(x)dx
and the inverse Laplace-Carson transform is defined as
(22)f(x)=Lc-1[Lc[f(x)]]=-12πi∫α-i∞α+i∞esx[-1[∫0sLc(t)+F(0)]]ds,
where F(s) is the Laplace transform of f(x).  Before we continue, we shall prove that the previous definition is indeed the inverse inverse Laplace-Carson. In fact, from the definition of inverse Laplace-Carson of a function f(x), we have that
(23)Lc(s)=Lc[f(x)](s)=∫0∞xe-xsf(x)dx=-dF(s)ds;
thus,
(24)∫0sLc(t)dt=-[F(s)-F(0)].
It follows that
(25)-12πi∫α-i∞α+i∞esx[-1[∫0sLc(t)dt+F(0)]]ds=-12πi∫α-i∞α+i∞esx[-[F(s)]]ds,Lc-1[Lc[f(x)]]=(-1)22πi∫α-i∞α+i∞esx[[F(s)]]ds=f(x).
Therefore, the inverse inverse Laplace-Carson is well defined.


## 5. Some Properties of Laplace-Carson Transform

In this part of the section, we consider some of the properties of the inverse Laplace-Carson that will enable us to find further transform pairs{f(x),Lc(s)} without having to compute and consider the following:
(26)(I)Lc[s+c]=Mn[e-cxf(x)](27)(II)Lc[f(ax)](s)=1aLc[sa](28)(III)∫α-i∞α+i∞esxLc(s)ds=xf(x)(29)(IV)Lc[af(x)+bg(x)](s)=[aLc(f(x))+bLc(g(x))](s)(30)(V)Lc[f(x)x](s)=L[f(x)](s)(31)(VI)Lc[f(x)*h(x)](s)=-[dF(s)dsG(s)+dG(s)dsF(s)](32)(VII)Lc[dnf(x)dxn](s)=-[dkf(0)dxknsn-1F(s)+sndF(s)ds-∑k=0n-2(n-k-1)sn-k-2dkf(0)dxk].
Let us verify the previous properties. We shall start with (I), by definition, we have the following:
(33)Lc[e-cxf(x)]=∫0∞[xe-cxe-sxf(x)]dx=∫0∞[xe-(c+s)xf(x)]dx=Lc[s+c]
and then the first property is verified.For (II) we have the following by definition:(34)Lc[f(ax)](s)=∫0∞[xe-xsf(ax)]dx=-dds[L[f(ax)](s)].
Now, using the property of the Laplace transform L[f(ax)](s)=(1/a)F(s/a), we can further obtain
(35)Lc[f(ax)](s)=-dds[1aF(sa)]=1a(-1)dds[F(sa)]=1aLc[sa]
and then, the property number (II) is verified.For number (III), we have the following: Letg(x)=xf(x); then
(36)∫α-i∞α+i∞esxLc(s)ds=∫α-i∞α+i∞esx[∫0∞e-xsxf(x)dx]ds=∫α-i∞α+i∞esx[∫0∞e-xsg(x)dx]ds.
By the theorem of inverse Laplace transform, we obtain
(37)∫α-i∞α+i∞esxMn(s)ds=g(x)=xf(x).
Number (IV) and (V) are obvious to be verified. For number (VI), we have the following by definition:
(38)Lc[f(x)*h(x)](s)=∫0∞[xe-sxf(x)*h(x)]=-dds[L(f(x)*h(x))(s)].
Now, using the property of Laplace transform of the convolution, we obtain the following:
(39)L(f(x)*h(x))(s)=F(s)·G(s)
and then, using the property of the derivative for the product of two functions, we obtain
(40)Lc[f(x)*h(x)](s)=-dds[F(s)·G(s)]=-[dF(s)dsG(s)+dG(s)dsF(s)].
And then, the property number (VI) is verified.For number (VII), by definition, we have the following:(41)Lc[dnf(x)dxn](s)=∫0∞[xe-sxdnf(x)dxn]dx=-dds[L(dnf(x)dxn)(s)].
Now, using the property of the Laplace transform,
(42)L(dnf(x)dxn)(s)=snF(s)-∑k=0n-1sn-k-1dkf(0)dxk.
Now, deriving the previous expression n-time, we obtain the following expression:
(43)-dds[snF(s)-∑k=0n-1sn-k-1dkf(0)dxk]=-[dkf(0)dxknsn-1F(s)+sndF(s)ds11111111-∑k=0n-2(n-k-1)sn-k-2dkf(0)dxk].
This completes the proof of number (VI). We shall now use some properties of the Laplace-Carson transform to solve the fractional groundwater flow equation. To achieve this, we shall start by assuming that the solution of the main equation can be separated as follows:
(44)Φ(r,t)=Φ1(t)Φ2(r);
then, the separated equations become
(45)DtαoCΦ1(t)+λ2Φ1(t)=0,(46)DrrΦ2(r)+1rDrΦ2(r)+λ2Φ2(r)=0,
where λ is the separation constant.The first equation (45) can be solve directly by applying on both sides the Laplace transform to obtain
(47)L(Φ1(t))=Θ1(s)=sα-1sα+λ2.
Using the inverse formula of Laplace transform of two-parameter Mittag-Leffler function, we get
(48)Φ1(t)=cEα,1(-STλ2tα),
where the Mittag-Leffler function is defined as follows:
(49)Eα,1(t)=∑n=0∞tαnΓ[αn+1].
To solve the second equation, we make use of the Laplace-Carson presented earlier, to obtain
(50)-Ds[s2Θ2(s)-sΦ2(0)-dΦ2(0)dr]+sΘ2(s)-Φ2(0)-DsΘ2(s)=0.
Deriving, we obtain the following ordinary differential equation:
(51)DsΘ2(s)[s2+λ2]=-sΘ2(s)
for which the exact solution is given as
(52)Θ2(s)=1s2+λ2.
Now applying the inverse Laplace transform operator on both sides in the previous equation, we obtain the following in terms of the Bessel function first kind [24]:
(53)Φ1(r)=J0(λr),
where the Bessel function first kind is defined as
(54)J0(r)=∑k=0∞(-1)kk!1Γ(k+1)(r2)2k.
Therefore, the solution of the fractional groundwater flow equation is given as
(55)Φ(r,t)=c∑n=0∞Eα,1(-STλn2tα)J0(λnr).
Making use of the initial and boundary conditions, we obtain the constant c to be
(56)c=Q4πT.
Then,
(57)Φ(r,t)=Q4πT∑n=0∞Eα,1(-STλn2tα)J0(λnr).
We shall present an alternative solution in the next section via the Boltzmann variable method. By using the boundary condition, we can determine the Eigen value of (55).


## 6. An Alternative Derivation of the Time-Fractional Groundwater Equation

In this section, we present an alternative approximate solution of the groundwater flow equation. A method frequently used to derive some kind of parabolic partial differential equations is the so-called Boltzmann transformation [5], defined for an arbitrary t0<t by  equation.
(58)u0=Sr24T(t-t0).
Let us consider now the following function:
(59)Φ(r,t)=c(t-t0)Eα,1[-u0]
with c an arbitrary constant.  If we assume that rb is the ratio of the borehole from which the water is taken out of the aquifer, then the total volume of the water withdrawn from the aquifer is given by
(60)Q0Δt0=4πcT.
Hence,
(61)Φ(r,t)=Q0Δt04πT(t-t0)Eα,1[-u0]
is the drawdown that will be observed at a distance, r, from the pumped borehole after the period Δt0.Now suppose that the a previous procedure is repeatedn times; that is, water is withdrawn for a short period of time, Δtk, at a consecutive times, tk+1=tk+Δtk, (k=0,1,…,n). Now, since the fractional groundwater flow equation is linear differential equation, the total drawdown at any time t>tn will be given by
(62)Φ(r,t)=14πT∑k=0nQkΔtk4πT(t-tk)Eα,1[-uk].
Therefore, if Δt→0, the definition of the defined integral can be invoked to write
(63)Φ(r,t)=14πT∫t0tQ(τ)dτ(t-τ)Eα,1[-Sr24T(t-τ)]dτ.
However, using the Boltzman variable, we arrive at the following expression:
(64)y=Sr24T(t-τ),Φ(r,t)=14πT∫t0∞Q(y)yEα,1[-y]dy.
The previous solution will be called the most generalized general solution of fractional groundwater equation. However, this solution can be simplified somewhat under certain conditions. A particularly important solution which arises when t0 is taken at zero and Q(t) is a constant independent of time, and then we arrive at
(65)Φ(r,t)=Q4πT∫u∞1yEα,1[-y]dy=Q4πTWα(u).
Here
(66)Wα(u)=∫u∞1yEα,1[-y]dy
will be called the generalized exponential integral.  It is worth pointing out that if alpha is equal to 1, we recover the exact analytical solution of the groundwater flow equation proposed by Theis. Alternative iterations method [8, 11, 20] can be used to derive approximate solutions of these problems.


## 7. Numerical Simulations

In this section, we investigate the behavior of the analytical solutions of the Theim fractional and the Theis fractional groundwater flow equation.  We compare the analytical solution with the experimental data from the pumping test obtained from the experimental site of the Institute for Groundwater Studies, the University of the Free State, Bloemfontein Campus, South Africa. We shall start with the simulation of Theim fractional groundwater flow equation. The analytical solution of the main problem was depicted in Figures1(a), 1(b), 1(c), and 1(d). It is worth to mention from the figures that the order of the derivative plays an important role in the simulation.Numerical simulation of the Theim-time-fractional order groundwater flow equation, ((a) showing the surface forα=0.95); ((b) α=0.75); ((c) α=0.65) and ((d) α=0.55).
(a)
(b)
(c)
(d)We shall present the numerical solution for the Theis-time-fractional groundwater flow equation for a fixed distance. The analytical solution of the main problem was depicted in Figures2(a), 2(b), 2(c), and 2(d).Numerical representation of (66) for different values of alpha.
(a)
(b)
(c)
(d)We shall present in Figure3 the comparison of the analytical solution of the Theis-time-fractional ground water flow equation, with the experimental data from the pumping test performed by the Institute for Groundwater Study on one of its borehole 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 Q=4.5 and monitoring the piezometric head for 350 minutes, at the distance of  r=31.54 m.Figure 3
Comparison of Theis-time-fractional groundwater flow equation with experimental data from real observation.A Chinese proverb says, I quote “An image is equivalent to ten thousand words.” With any reservation we can conclude that the fractional order derivative plays an important role while modeling real-world problem into mathematical equation. 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 borehole settled on the test site of the University of the Free State.


## 8. Conclusion

We have generalized to Theim and Theis groundwater flow equation to the concept of fractional order derivatives.  The resulting equations were solved analytically via several techniques, including the Laplace transform method, the Fourier transform method, the Laplace-Carson transform method, and the Boltzmann variable method. The numerical simulations show that the fractional order derivative plays an important role in the simulation process. In addition, we compare the analytical solution with experimental data to access the accuracy of the fractional groundwater model. The analytical solution was in perfect agreement with experimental data.



---

*Source: 543026-2013-10-27.xml*