The Use of Fractional Order Derivative to Predict the Groundwater Flow

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 rate Q and monitoring the piezometric head for 350 minutes.

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Γ()𝑑𝑛−𝛼𝑑𝑥^𝑛𝑛×^𝑥₀(𝑥−𝑡)^{𝑛−𝛼−1}𝑓(𝑡)𝑑𝑡,𝑛−1≤𝛼≤𝑛. (1) For the case of the Caputo fractional order derivative, we have the following definition: 𝐷^𝛼_𝑥𝑓(𝑥=1)Γ()×𝑛−𝛼_𝑥0(𝑥−𝑡)^{𝑛−𝛼−1}𝑑^𝑛𝑑𝑥^𝑛𝑓(𝑡)𝑑𝑡,𝑛−1≤𝛼≤𝑛. (2) 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]) 𝐷^𝛼_𝑥𝑓(𝑥)=1Γ()𝑑𝑛−𝛼𝑑𝑥^𝑛𝑛×^𝑥₀(𝑥−𝑡)^{𝑛−𝛼−1}𝑓(𝑡)−𝑓(0)𝑑𝑡. (3) 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]: 𝐷^𝛼_𝑟𝑟Φ(𝑟)+1𝑟Φ(𝑟)=0,1<𝛼≤2. (4) Subject to the initial condition, 𝑄=2𝜋𝑇𝐷_𝑟(Φ(𝑟_𝑏)).

We will make use of the Laplace transform to derive analytical solution of (4). Thus, multiplying on both sides of (4) by 𝑟 and secondly applying the Laplace transform, we obtain the following expression: 𝑑[𝐿(Φ)(𝑠)]+𝛼𝑑𝑠𝑠+1𝑠^𝛼(𝐿(Φ𝑠=)())_𝑙𝑚=2𝑑_𝑚(𝑚−1)𝑠^{𝑚−2−𝛼}, (5) where 𝑑_𝑚=𝐷₀^𝛼−𝑚+Φ(0⁺)(𝑚=2,…,𝑙). Now, one can derive the solution of the ordinary order differential equation with respect to the Laplace transform of Θ(𝑠)=𝐿(Φ(𝑟)): Θ(𝑠)=𝑠^−𝛼−𝑠exp^1−𝛼×𝑎1−𝛼₁+^𝑙_𝑚=2𝑑_𝑚(𝑚−1)𝑠^𝑚−2−𝑠exp^1−𝛼,1−𝛼𝑑𝑠 (6) with 𝑎₁ 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: Θ(𝑠)=𝑐Θ₁(𝑠)+^𝑙_𝑚=2𝑑_𝑚(𝑚−1)Θ^∗_𝑚(𝑠) (7) with, of course, Θ₁(𝑠)=𝑠^−𝛼−𝑠exp^1−𝛼,Θ1−𝛼_∗𝑚(𝑠)=𝑠^−𝛼𝑠exp^1−𝛼𝛼−1^∞_𝑗=011−𝛼^𝑗𝑠^{(1−𝛼)𝑗+𝑚−1}).(1−𝛼𝑗+𝑚−1𝑗! (8) Now, applying the inverse Laplace transform on Θ₁(𝑠) and using the fact that 𝑠^{−[𝛼+(𝛼−1)𝑗]}𝑟=𝐿^{𝛼+(𝛼−1)𝑗−1}Γ𝛼+(𝛼−1)𝑗 (9) we obtain Φ₁(𝑟)=𝑟^𝛼−1𝑜Ψ₁)∣𝑥(𝛼,𝛼−1^𝛼−1𝛼−1 (10) with 𝑜Ψ₁[⋅] the generalized Wright function [16] for 𝑝=1 and 𝑞=2. We next expand the exponential function exp[−𝑠^1−𝛼/(1−𝛼)] in power series; multiplying the resulting two series; in addition of this if we consider the number 𝑏_𝑘(𝛼,𝑚) defined for 𝛼>0, 𝑚=2,…𝑙, (𝛼≠(𝑝+𝑚−1)/𝑝,𝑝∉ℕ) and 𝑘∈ℕ₀, 𝑏_𝑘(𝛼,𝑚)=^𝑙_{𝑝,𝑗=0,…𝑘,𝑝+𝑗=𝑘}(−1)^𝑞().𝑝!𝑗!1−𝛼𝑞+𝑚−1 (11) The previous family of number possesses satisfies the following recursive formula: 𝑏_𝑘(𝛼,𝑚)𝑏_𝑘+1()=𝛼,𝑚𝛼−𝑚𝛼−1+𝑘, (12) which produces the explicit expression for 𝑏_𝑘(𝛼,𝑚) in the form of 𝑏_𝑘()=Γ)/(𝛼,𝑚[(𝛼−𝑚𝛼−1)](𝑚−1)Γ[((𝛼−𝑚)/(𝛼−1))+𝑘],𝑘∈ℕ₀. (13) Now, having the previous expression on hand, we can derive that Θ^∗_𝑚(𝑠)=𝑠^{𝑚−𝛼−1∞}_𝑗=011−𝛼^𝑗𝑠^{(1−𝛼)𝑝}×𝑝!_∞𝑝=011−𝛼^𝑝(−1)^𝑝)𝑠(1−𝛼𝑝+𝑚−1^{(1−𝛼)𝑗}=𝑝!_∞𝑘=0𝑏_𝑘()1𝛼,𝑚1−𝛼^𝑘𝑠^{(1−𝛼)𝑘+𝑚−𝛼−1}(𝑚=2,…,𝑙). (14) However, remembering (9) with 𝛽=(𝛼−1)𝑘+𝛼+1−𝑚, we can further derive the following expression for Φ^∗_𝑚(𝑟) as Φ^∗_𝑚(𝑟)=^∞_𝑘=0𝑏_𝑘(1𝛼,𝑚)1−𝛼^𝑘×Γ(𝑘+1)Γ[()𝑘]𝑥𝛼+1−𝑚+𝛼−1𝑘!^{(𝛼−1)𝑘+𝛼−𝑚} (15) or in the simplified version we have Φ^∗_𝑚(𝑟)=Γ[(𝛼−𝑚)/(𝛼−1)]()Φ𝑚−1^𝑚(𝑟), (16) where Φ_𝑚(𝑟)=𝑟₁^𝛼−𝑚Ψ₂×()(),1,1𝛼+1−𝑚,𝛼−1𝛼−𝑚∣𝑟𝛼−1,1^𝛼−1.𝛼−1 (17) It follows that the solution of the fractional-Thiem groundwater flow equation is in the form of Φ(𝑟)=𝑎₁𝑟^𝛼−1𝑜Ψ₁)∣𝑥(𝛼,𝛼−1^𝛼−1𝛼−1+𝑎₂^𝑙_𝑚=2𝑏_𝑚()Γ)/(𝑚−1[(𝛼−𝑚𝛼−1)]()𝑟𝑚−1^𝛼−𝑚×₁Ψ₂()(),1,1𝛼+1−𝑚,𝛼−1𝛼−𝑚∣𝑟𝛼−1,1^𝛼−1.𝛼−1 (18) In our case 𝑙=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: 𝑆𝐷_𝑡Φ(𝑟,𝑡)=𝑇𝐷_𝑟𝑟Φ()+1𝑟,𝑡𝑟𝐷_𝑟Φ(𝑟,𝑡). (19) 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 𝑆𝐷^𝛼_𝑡Φ(𝑟,𝑡)=𝑇𝐷_𝑟𝑟Φ()+1𝑟,𝑡𝑟𝐷_𝑟Φ(𝑟,𝑡),0<𝛼≤1 (20) with initial condition Φ(𝑟,0)=0 and boundary condition lim_𝑟→∞Φ(𝑟,𝑡)=0, 𝑄=2𝜋𝑇𝜕𝑟Φ(𝑟_𝑏,𝑡), here 𝑇 is the transmissivity of the aquifer, 𝑟_𝑏 is the ratio of the borehole, and 𝑄 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. Let 𝑓(𝑥) be a continuous function over an open interval (0,∞) such that its Laplace transform is 𝑛 time differentiable; then the Laplace-Carson transform 𝑓 is defined as follows: 𝐿_𝑐(𝑠)=𝐿_𝑐𝑓(𝑥)(𝑠)=^∞₀𝑥𝑒^−𝑥𝑠𝑓(𝑥)𝑑𝑥 (21) and the inverse Laplace-Carson transform is defined as 𝑓(𝑥)=𝐿_𝑐⁻¹𝐿_𝑐𝑓(𝑥=)−12𝜋𝑖^𝛼+𝑖∞_{𝛼−𝑖∞}𝑒^𝑠𝑥−1^𝑠₀𝐿_𝑐(𝑡)+𝐹(0)𝑑𝑠, (22) where 𝐹(𝑠) is the Laplace transform of 𝑓(𝑥). 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 𝑓(𝑥), we have that 𝐿_𝑐(𝑠)=𝐿_𝑐𝑓(𝑥)(𝑠)=^∞₀𝑥𝑒^−𝑥𝑠𝑓(𝑥)(𝑠)𝑑𝑥=−𝑑𝐹;𝑑𝑠 (23) thus, ^𝑠₀𝐿_𝑐(𝑡)𝑑𝑡=−[𝐹(𝑠)−𝐹(0)]. (24) It follows that −12𝜋𝑖^𝛼+𝑖∞_{𝛼−𝑖∞}𝑒^𝑠𝑥−1^𝑠₀𝐿_𝑐(𝑡)(0=𝑑𝑡+𝐹)𝑑𝑠−12𝜋𝑖^𝛼+𝑖∞_{𝛼−𝑖∞}𝑒^𝑠𝑥[−[𝐹(𝑠𝐿)]]𝑑𝑠,_𝑐−1𝐿_𝑐𝑓(𝑥=())−1²2𝜋𝑖^𝛼+𝑖∞_{𝛼−𝑖∞}𝑒^𝑠𝑥[[𝐹(𝑠)]]𝑑𝑠=𝑓(𝑥). (25) 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 {𝑓(𝑥),𝐿_𝑐(𝑠)} without having to compute and consider the following: (I)𝐿_𝑐[𝑠+𝑐]=𝑀_𝑛𝑒^−𝑐𝑥𝑓(𝑥()𝐿)II_𝑐𝑓(𝑠)=1𝑎𝑥)(𝑎𝐿_𝑐𝑠𝑎(III)^𝛼+𝑖∞_{𝛼−𝑖∞}𝑒^𝑠𝑥𝐿_𝑐(𝑠)(𝑥)()𝐿𝑑𝑠=𝑥𝑓IV_𝑐(𝑥)(𝑥𝑠)=𝑎𝑓+𝑏𝑔)(𝑎𝐿_𝑐𝑓(𝑥)+𝑏𝐿_𝑐𝑔(𝑥)𝑠)(V)𝐿_𝑐(𝑓(𝑥)𝑥𝑠)𝑓(𝑥𝑠)()𝐿(=𝐿)(VI_𝑐𝑓(𝑥)(𝑥𝑠)(𝑠)∗ℎ)(=−𝑑𝐹𝐺(𝑠)+(𝑠)𝑑𝑠𝑑𝐺𝐹(𝑠()𝐿𝑑𝑠)VII_𝑐𝑑^𝑛𝑓(𝑥)𝑑𝑥^𝑛(𝑠)=−𝑛𝑠^𝑛−1𝐹(𝑠)+𝑠^𝑛𝑑𝐹(𝑠)−𝑑𝑠_{𝑛−2𝑘=0}(𝑛−𝑘−1)𝑠^{𝑛−𝑘−2}𝑑^𝑘𝑓(0)𝑑𝑥^𝑘. (26)(27)(28)(29)(30)(31)(32) Let us verify the previous properties. We shall start with (I), by definition, we have the following: 𝐿_𝑐𝑒^−𝑐𝑥𝑓(𝑥)=^∞₀𝑥𝑒^−𝑐𝑥𝑒^−𝑠𝑥𝑓(𝑥=)𝑑𝑥_∞0𝑥𝑒^{−(𝑐+𝑠)𝑥}𝑓(𝑥)𝑑𝑥=𝐿_𝑐[𝑠+𝑐] (33) and then the first property is verified.

For (II) we have the following by definition: 𝐿_𝑐𝑓(𝑠)=𝑎𝑥)(_∞0𝑥𝑒^−𝑥𝑠𝑓(𝑑𝑎𝑥)𝑑𝑥=−𝐿𝑓(𝑠.𝑑𝑠𝑎𝑥)() (34) Now, using the property of the Laplace transform 𝐿[𝑓(𝑎𝑥)](𝑠)=(1/𝑎)𝐹(𝑠/𝑎), we can further obtain 𝐿_𝑐𝑓(𝑠)𝑑𝑎𝑥)(=−1𝑑𝑠𝑎𝐹𝑠𝑎=1𝑎()𝑑−1𝐹𝑠𝑑𝑠𝑎=1𝑎𝐿_𝑐𝑠𝑎 (35) and then, the property number (II) is verified.

For number (III), we have the following: Let 𝑔(𝑥)=𝑥𝑓(𝑥); then ^𝛼+𝑖∞_{𝛼−𝑖∞}𝑒^𝑠𝑥𝐿_𝑐(𝑠)𝑑𝑠=^𝛼+𝑖∞_{𝛼−𝑖∞}𝑒^𝑠𝑥∞₀𝑒^−𝑥𝑠(𝑥)=𝑥𝑓𝑑𝑥𝑑𝑠_{𝛼+𝑖∞𝛼−𝑖∞}𝑒^𝑠𝑥∞₀𝑒^−𝑥𝑠𝑔(𝑥)𝑑𝑥𝑑𝑠. (36) By the theorem of inverse Laplace transform, we obtain ^𝛼+𝑖∞_{𝛼−𝑖∞}𝑒^𝑠𝑥𝑀_𝑛(𝑠)𝑑𝑠=𝑔(𝑥)=𝑥𝑓(𝑥). (37) Number (IV) and (V) are obvious to be verified. For number (VI), we have the following by definition: 𝐿_𝑐𝑓(𝑥)(𝑥𝑠)=∗ℎ)(_∞0𝑥𝑒^−𝑠𝑥𝑓(𝑥)(𝑥𝑑∗ℎ)=−𝐿𝑓(𝑥)(𝑥𝑠.𝑑𝑠∗ℎ)() (38) Now, using the property of Laplace transform of the convolution, we obtain the following: 𝐿𝑓(𝑥)∗ℎ(𝑥)(𝑠)=𝐹(𝑠)⋅𝐺(𝑠) (39) and then, using the property of the derivative for the product of two functions, we obtain 𝐿_𝑐𝑓(𝑥)(𝑥𝑠)𝑑∗ℎ)(=−[𝐹(𝑠)(𝑠(𝑠)𝑑𝑠⋅𝐺)]=−𝑑𝐹𝐺(𝑠)+(𝑠)𝑑𝑠𝑑𝐺𝐹(𝑠.𝑑𝑠) (40) And then, the property number (VI) is verified.

For number (VII), by definition, we have the following: 𝐿_𝑐𝑑^𝑛𝑓(𝑥)𝑑𝑥^𝑛(𝑠)=^∞₀𝑥𝑒^−𝑠𝑥𝑑^𝑛𝑓(𝑥)𝑑𝑥^𝑛𝑑𝑑𝑥=−𝐿𝑑𝑑𝑠^𝑛𝑓(𝑥)𝑑𝑥^𝑛(𝑠). (41) Now, using the property of the Laplace transform, 𝐿𝑑^𝑛𝑓(𝑥)𝑑𝑥^𝑛(𝑠)=𝑠^𝑛𝐹(𝑠)−^𝑛−1_𝑘=0𝑠^{𝑛−𝑘−1}𝑑^𝑘𝑓(0)𝑑𝑥^𝑘. (42) Now, deriving the previous expression n-time, we obtain the following expression: −𝑑𝑠𝑑𝑠^𝑛𝐹(𝑠)−^𝑛−1_𝑘=0𝑠^{𝑛−𝑘−1}𝑑^𝑘𝑓(0)𝑑𝑥^𝑘=−𝑛𝑠^𝑛−1𝐹(𝑠)+𝑠^𝑛𝑑𝐹(𝑠)−𝑑𝑠_{𝑛−2𝑘=0}(𝑛−𝑘−1)𝑠^{𝑛−𝑘−2}𝑑^𝑘𝑓(0)𝑑𝑥^𝑘. (43) 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) then, the separated equations become ^𝐶_𝑜𝐷^𝛼_𝑡Φ₁(𝑡)+𝜆²Φ₁(𝑡)𝐷=0,_𝑟𝑟Φ₂(𝑟)+1𝑟𝐷_𝑟Φ₂(𝑟)+𝜆²Φ₂(𝑟)=0, (45)(46) where 𝜆 is the separation constant.

The first equation (45) can be solve directly by applying on both sides the Laplace transform to obtain 𝐿Φ₁(𝑡)=Θ₁(𝑠)=𝑠^𝛼−1𝑠^𝛼+𝜆². (47) Using the inverse formula of Laplace transform of two-parameter Mittag-Leffler function, we get Φ₁(𝑡)=𝑐𝐸_𝛼,1−𝑆𝑇𝜆²𝑡^𝛼, (48) where the Mittag-Leffler function is defined as follows: 𝐸_𝛼,1(𝑡)=^∞_𝑛=0𝑡^𝛼𝑛Γ[].𝛼𝑛+1 (49) To solve the second equation, we make use of the Laplace-Carson presented earlier, to obtain −𝐷_𝑠𝑠²Θ₂(𝑠)−𝑠Φ₂(0)−𝑑Φ₂(0)𝑑𝑟+𝑠Θ₂(𝑠)−Φ₂(0)−𝐷_𝑠Θ₂(𝑠)=0. (50) Deriving, we obtain the following ordinary differential equation: 𝐷_𝑠Θ₂(𝑠)𝑠²+𝜆²=−𝑠Θ₂(𝑠) (51) for which the exact solution is given as Θ₂(𝑠)=1𝑠²+𝜆². (52) 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) where the Bessel function first kind is defined as 𝐽₀(𝑟)=^∞_𝑘=0(−1)^𝑘1𝑘!Γ()𝑟𝑘+12^2𝑘. (54) Therefore, the solution of the fractional groundwater flow equation is given as Φ(𝑟,𝑡)=𝑐^∞_𝑛=0𝐸_𝛼,1−𝑆𝑇𝜆²_𝑛𝑡^𝛼𝐽₀𝜆_𝑛𝑟. (55) Making use of the initial and boundary conditions, we obtain the constant 𝑐 to be 𝑄𝑐=.4𝜋𝑇 (56) Then, Φ()=𝑄𝑟,𝑡4𝜋𝑇^∞_𝑛=0𝐸_𝛼,1−𝑆𝑇𝜆²_𝑛𝑡^𝛼𝐽₀𝜆_𝑛𝑟. (57) 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 𝑡₀<𝑡 by equation. 𝑢₀=𝑆𝑟²4𝑇𝑡−𝑡₀. (58) Let us consider now the following function: Φ()=𝑐𝑟,𝑡𝑡−𝑡₀𝐸_𝛼,1−𝑢₀ (59) with 𝑐 an arbitrary constant. If we assume that 𝑟_𝑏 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 𝑄₀Δ𝑡₀=4𝜋𝑐𝑇. (60) Hence, Φ()=𝑄𝑟,𝑡⁰Δ𝑡₀4𝜋𝑇𝑡−𝑡₀𝐸_𝛼,1−𝑢₀ (61) is the drawdown that will be observed at a distance, 𝑟, from the pumped borehole after the period Δ𝑡₀.

Now suppose that the a previous procedure is repeated 𝑛 times; that is, water is withdrawn for a short period of time, Δ𝑡_𝑘, at a consecutive times, 𝑡_𝑘+1=𝑡_𝑘+Δ𝑡_𝑘, (𝑘=0,1,…,𝑛). Now, since the fractional groundwater flow equation is linear differential equation, the total drawdown at any time 𝑡>𝑡_𝑛 will be given by Φ()=1𝑟,𝑡4𝜋𝑇^𝑛_𝑘=0𝑄_𝑘Δ𝑡_𝑘4𝜋𝑇𝑡−𝑡_𝑘𝐸_𝛼,1−𝑢_𝑘. (62) Therefore, if Δ𝑡→0, the definition of the defined integral can be invoked to write Φ()=1𝑟,𝑡4𝜋𝑇^𝑡_𝑡0𝑄(𝜏)𝑑𝜏()𝐸𝑡−𝜏^𝛼,1−𝑆𝑟²4𝑇(𝑡−𝜏)𝑑𝜏. (63) However, using the Boltzman variable, we arrive at the following expression: 𝑦=𝑆𝑟²(),Φ()=14𝑇𝑡−𝜏𝑟,𝑡4𝜋𝑇^∞_𝑡0𝑄𝑦𝑦𝐸_𝛼,1−𝑦𝑑𝑦. (64) 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 𝑡₀ is taken at zero and 𝑄(𝑡) is a constant independent of time, and then we arrive at Φ()=𝑄𝑟,𝑡4𝜋𝑇^∞_𝑢1𝑦𝐸_𝛼,1𝑄−𝑦𝑑𝑦=𝑊4𝜋𝑇^𝛼(𝑢). (65) Here 𝑊_𝛼(𝑢)=^∞_𝑢1𝑦𝐸_𝛼,1−𝑦𝑑𝑦 (66) 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 Figures 1(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.

(a)

(b)

(c)

(d)

Figure 1 

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).

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 Figures 2(a), 2(b), 2(c), and 2(d).

(a)

(b)

(c)

(d)

Figure 2 

Numerical representation of (66) for different values of alpha.

We shall present in Figure 3 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 𝑄=4.5 and monitoring the piezometric head for 350 minutes, at the distance of 𝑟=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.