# A Modified Groundwater Flow Model Using the Space Time Riemann-Liouville Fractional Derivatives Approximation

Authors: Abdon Atangana; S. C. Oukouomi Noutchie
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

The notion of uncertainty in groundwater hydrology is of great importance as it is known to result in misleading output when neglected or not properly accounted for. In this paper we examine this effect in groundwater flow models. To achieve this, we first introduce the uncertainties functionsu as function of time and space. The function u accounts for the lack of knowledge or variability of the geological formations in which flow occur (aquifer) in time and space. We next make use of Riemann-Liouville fractional derivatives that were introduced  by Kobelev and Romano in 2000 and its approximation to modify the standard version of groundwater flow equation. Some properties of the modified Riemann-Liouville fractional derivative approximation are presented. The classical model for groundwater flow, in the case of density-independent flow in a uniform homogeneous aquifer is reformulated by replacing the classical derivative by the Riemann-Liouville fractional derivatives approximations. The modified equation is solved via the technique of green function and the variational iteration method.


---

## Body

## 1. Introduction

This paper investigates the effects of uncertainty on the predictive accuracy of flow through porous media; it is commonly believed that the problem that occurs in groundwater models is the suitable geometry in which flow occurs on one hand and the deviation of theoretical expected values from observations on the other hand. Therefore, it is important to notice that miniscule effects observed always require the most new modifications of ideas. Scientists in the field of hydrogeology in particular are used to deal with doubt and uncertainty, because it is impossible to understand or to model the phenomena that occur in aquifers exactly. All historical and current theoretical knowledge in groundwater investigations are uncertain and doubtful. This experience with doubt and uncertainty is important. We believe that it is of great importance, and one that extends beyond the theories which are used to interpret the phenomena that take place in aquifers. Doubt is clearly a value that must be analytically included in groundwater flow models. Uncertainty in groundwater hydrology originates from different sources. Neglecting uncertainty in groundwater assessments can lead to incorrect results and misleading output. Generally there are various sources of uncertainty in model outputs, for example, uncertainty associated with lack of knowledge or accuracy of the model inputs as well as the structural uncertainty related to the mathematical interpretation of the model. The assessment and presentation of the effects of uncertainty are now widely recognized as important parts of analyses for complex systems [1–6]. At the simplest level such analyses can be viewed as the study of functions. In order to include explicitly the possible effect of the uncertainties into mathematical models, we introduced in this paper the uncertainties in groundwater models as a function of time and space.Consider(1)
u
=
u
(
x
_
,
t
)
.


## 2. Modification of Groundwater Flow Equation

To be clear, the modification of the classical model for groundwater flow in the case of density independent flow in the uniform and homogeneous aquifer is considered in this paper.To modify this, we make use of Riemann-Liouville fractional derivatives that were introduced in [7] and attempted by many others, see for example, [8]. These derivatives are defined as
(2)
D
1
+
ε
t
f
=
D
+
,
t
μ
t
f
=
(
d
d
t
)
n
∫
0
t
[
f
(
τ
)
Γ
(
n
-
μ
τ
(
τ
)
(
t
-
τ
)
μ
τ
-
n
+
1
)
]
d
τ
D
1
+
ε
x
f
=
D
+
,
x
μ
x
f
=
(
d
d
x
)
n
∫
0
x
[
f
(
τ
)
Γ
(
n
-
μ
τ
(
τ
)
(
x
-
τ
)
μ
τ
-
n
+
1
)
]
d
τ
.
Here, Γ is the Euler gamma function; n
=
⌈
μ
⌉+1, where ⌈
μ
⌉ is the integer part of μ for μ
≥
0 that is n
-
1
≤
μ
<
n and n
=
0 for μ
<
n. Following equation (2) we have that μ
t
=
1
+
ε
t and μ
x
=
1
+
ε
x. The integral operator defined above for fractional exponents μ
x and μ
t depending on coordinates and time can be expressed in terms of ordinary derivative and integral [7] for |
ε
|
≪
1. For this matter, generalized Riemann-Liouville fractional derivatives satisfy the approximate relations.Consider(3)
D
1
+
u
t
f
≅
(
1
+
u
t
)
∂
f
∂
t
+
∂
u
t
∂
t
f
D
1
+
u
x
f
≅
(
1
+
u
x
)
∂
f
∂
x
+
∂
u
x
∂
x
f
.
The above relations make it possible to describe the flow system, including the effect of uncertainties on the behaviour of physical systems, by means of partial differential and integral equations.Let us examine some properties of the above derivative operator [9].(i)
Addition.
Ifu
x, f
(
x
) and g
(
x
) are differentiable in the opened intervalI, then
(4)
D
1
+
u
x
[
f
(
x
)
+
g
(
x
)
]
≅
D
1
+
u
x
[
f
(
x
)
]
+
D
1
+
u
x
[
g
(
x
)
]
D
1
+
u
x
[
f
(
x
)
+
g
(
x
)
]
≅
(
1
+
u
x
)
∂
[
f
(
x
)
+
g
(
x
)
]
∂
x
+
∂
u
x
∂
x
[
f
(
x
)
+
g
(
x
)
]
(
1
+
u
x
)
∂
[
f
(
x
)
]
∂
x
+
∂
u
x
∂
x
[
f
(
x
)
]
+
(
1
+
u
x
)
∂
[
g
(
x
)
]
∂
x
+
∂
u
x
∂
x
[
g
(
x
)
]
≅
D
1
+
u
x
[
f
(
x
)
]
+
D
1
+
u
x
[
f
(
x
)
]
.
(ii)
Division.
Ifu
x and 1
/
f
(
x
) are differentiable on the opened intervalI, then
(5)
D
1
+
u
x
[
1
f
(
x
)
]
≅
[
-
(
1
+
u
x
)
f
′
(
x
)
+
u
x
′
f
(
x
)
]
f
2
(
x
)
=
-
f
′
(
x
)
f
2
(
x
)
-
u
x
f
′
(
x
)
f
2
(
x
)
+
u
x
′
f
(
x
)
f
2
(
x
)
.
(iii)
Multiplication.
Ifu
x, f
(
x
) and g
(
x
) are differentiable in the opened intervalI, then
(6)
D
1
+
u
x
[
f
(
x
)
·
g
(
x
)
]
≅
g
(
x
)
f
′
(
x
)
+
f
(
x
)
g
′
(
x
)
+
(
g
f
′
+
f
g
′
)
(
x
)
u
x
+
u
x
′
(
f
(
x
)
g
(
x
)
)
.
(iv)
Power.
Ifu
x, and f
(
x
) are differentiable in the opened intervalI, then
(7)
D
1
+
u
x
[
(
f
(
x
)
)
n
]
≅
n
f
′
f
n
-
1
+
u
x
n
f
′
f
n
-
1
+
u
x
′
f
n
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
n
≥
1
.
(v)
Ifu
x, and f
(
x
) are two times differentiable in the opened intervalI, then
(8)
D
1
+
u
x
[
D
1
+
u
x
[
f
(
x
)
]
]
≅
(
1
+
u
x
)
[
(
1
+
u
x
)
∂
2
f
∂
x
2
+
3
∂
f
∂
x
∂
u
x
∂
x
+
∂
2
u
x
∂
x
2
f
]
+
∂
u
x
∂
x
f
.It is important to observe that if u
x
=
0, we recover the properties of normal derivatives. Recent investigations suggest that the flow is influenced by the geometry of the bedding parallel factures. An attempt to circumvent this problem, Barker introduced a model in which the geometry of the aquifer is regarded as a fractal [10]. In the same direction, the authors in [11] introduced the concept of noninteger 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 as results of their investigation, they found that there is a close relationship between the fractal and the fractional. Therefore, to include the fractal dimension into the mathematical formulation of the modified groundwater flow equation we next introduce the constant fractal dimension α. The classical model for groundwater flow, in the case of density-independent flow in a uniform homogeneous aquifer, can then be reformulated as follows:
(9)
D
1
+
u
t
(
r
,
t
)
Φ
(
r
,
t
)
=
K
S
0
D
1
+
u
t
(
r
,
t
)
[
r
D
1
+
u
t
(
r
,
t
)
Φ
(
r
,
t
)
]
+
f
(
r
,
t
)
,
where K is the hydraulic conductivity of the aquifer, S
0 is the specific storativity of the aquifer, f is the strength of any sources or sink, here it will be neglected, and finally Φ
(
r
,
t
) is the piezometric head.In order to meet the physical and mathematical requirements we impose the uncertainties function to be a positive function such that(10)
0
<
u
x
<
1
.
Equation (9) makes it possible to describe the flow through the geological formation, and the effect of uncertainties on the behaviour of physical systems, by means of partial differential and integral equations. However there is no analytical solution for this equation, in fact the analytical solution is very difficult to determine. Therefore we need the following approximation to simplify (9):
(11)
D
1
+
u
t
(
r
,
t
)
[
r
D
1
+
u
t
(
r
,
t
)
Φ
(
r
,
t
)
]
≅
(
1
+
u
r
)
[
r
(
1
+
u
r
)
∂
2
Φ
(
r
,
t
)
∂
r
2
+
r
∂
2
u
r
∂
r
2
Φ
(
r
,
t
)
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
∂
u
r
∂
r
Φ
(
r
,
t
)
+
(
1
+
u
r
)
∂
Φ
(
r
,
t
)
∂
r
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
3
r
∂
u
r
∂
r
∂
Φ
(
r
,
t
)
∂
r
∂
2
Φ
(
r
,
t
)
∂
r
2
+
r
∂
2
u
r
∂
r
2
Φ
(
r
,
t
)
+
∂
u
r
∂
r
]
+
r
(
∂
u
r
∂
r
)
2
Φ
(
r
,
t
)
Making use of (9), (3), and (11) we obtain the following equation:
(12)
(
1
+
u
t
)
∂
Φ
(
r
,
t
)
∂
t
+
∂
u
t
Φ
(
r
,
t
)
∂
t
≅
K
r
S
0
{
(
1
+
u
r
)
[
r
(
1
+
u
r
)
∂
2
Φ
(
r
,
t
)
∂
r
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
r
∂
2
u
r
∂
r
2
Φ
(
r
,
t
)
+
∂
u
r
∂
r
Φ
(
r
,
t
)
h
h
h
h
h
h
h
h
h
h
h
h
h
h
∂
2
Φ
(
r
,
t
)
∂
r
2
+
(
1
+
u
r
)
∂
Φ
(
r
,
t
)
∂
r
+
3
r
∂
u
r
∂
r
∂
Φ
(
r
,
t
)
∂
r
]
h
h
h
h
h
h
h
+
r
(
∂
u
r
∂
r
)
2
Φ
(
r
,
t
)
}
.
Since uncertainties additions to unit are small, the right- and left-hand sides of (12) can be divided by (1
+
u
t) to obtain the following approximate equation:(13a)
∂
Φ
(
r
,
t
)
∂
t
=
K
S
0
{
1
r
∂
Φ
(
r
,
t
)
∂
r
+
∂
2
Φ
(
r
,
t
)
∂
r
2
h
h
h
h
+
(
2
u
r
-
u
t
)
∂
2
Φ
(
r
,
t
)
∂
r
2
h
h
h
h
+
(
1
+
u
r
-
u
t
)
∂
2
u
r
∂
r
2
Φ
(
r
,
t
)
h
h
h
h
+
u
r
-
u
t
r
∂
u
r
∂
r
Φ
(
r
,
t
)
h
h
h
h
+
2
u
r
-
u
t
r
∂
Φ
(
r
,
t
)
∂
r
h
h
h
h
+
3
(
1
+
u
r
-
u
t
)
∂
u
r
∂
r
∂
Φ
(
r
,
t
)
∂
r
h
h
h
h
+
(
1
+
u
r
-
u
t
)
(
∂
u
r
∂
r
)
2
Φ
(
r
,
t
)
}
-
u
t
∂
u
t
∂
t
Φ
(
r
,
t
)
,
where u
r
2 term is omitted because it is significantly very small and since we are dealing with approximation here, it needs not to be considered in this case. For simplicity (13a) can be reformulated as
(13b)
∂
Φ
(
r
,
t
)
∂
t
=
K
S
0
[
1
r
∂
Φ
(
r
,
t
)
∂
r
+
∂
2
Φ
(
r
,
t
)
∂
r
2
]
+
F
(
u
r
,
u
t
,
Φ
(
r
,
t
)
)
.Here the additional term can be roughly approximate to(14)
F
(
u
r
,
u
t
,
Φ
(
r
,
t
)
)
=
[
(
2
u
r
-
u
t
)
∂
2
Φ
(
r
,
t
)
∂
r
2
h
+
3
(
1
+
u
r
-
u
t
)
∂
u
r
∂
r
∂
Φ
(
r
,
t
)
∂
r
h
-
u
t
(
∂
2
u
r
∂
r
2
+
1
r
∂
u
r
∂
r
)
Φ
(
r
,
t
)
h
∂
2
Φ
(
r
,
t
)
∂
r
2
+
u
r
(
∂
u
t
∂
r
+
∂
u
t
∂
t
)
Φ
(
r
,
t
)
]
and Φ
(
r
,
t
) satisfies the equation of classical model for groundwater flow in the case of density independent flow in the uniform and homogeneous aquifer. It is important to observe that the modified equations (13a) and (13b) differ from the standard form of groundwater flow equation in three properties:There is a new operator that takes into account the variation in piezometric head and uncertainties function given below as(15)
ℶ
(
u
r
,
u
t
,
Φ
)
=
(
2
u
r
-
u
t
)
∂
2
Φ
(
r
,
t
)
∂
r
2
+
3
(
1
+
u
r
-
u
t
)
∂
u
r
∂
r
∂
Φ
(
r
,
t
)
∂
r
.
Second, the “force”
(16)
F
=
∂
2
u
r
∂
r
2
+
1
r
∂
u
r
∂
r
appears due to the coordinate dependence of uncertainty function. And finally, there is a derivative-free term that depends only on the uncertainties time function
(17)
B
=
∂
u
t
∂
r
+
∂
u
t
∂
t
and is proportional to the piezometric head Φ and characterises, depending on the coefficient sign, the retardation or enhancement of the flow through the porous media.It is important to point out that those terms in (13a) that involve fractional additions, F and B, to the time and space dimensions are small. It follows that this equation can be solved approximately by changing the function Φby Φ
0, which satisfies the standard version of the groundwater flow equation which is the left-side of (13a) and (13b), in terms concerning u. Now let us suppose that such change is made in the expression (14). Equations (13a) and (13b) become
(18)
∂
Φ
(
r
,
t
)
∂
t
=
K
S
0
[
1
r
∂
Φ
(
r
,
t
)
∂
r
+
∂
2
Φ
(
r
,
t
)
∂
r
2
]
+
[
(
2
u
r
-
u
t
)
∂
2
Φ
(
r
,
t
)
∂
r
2
+
3
(
1
+
u
r
-
u
t
)
h
h
h
h
h
h
h
×
∂
u
r
∂
r
∂
Φ
(
r
,
t
)
∂
r
-
u
t
(
∂
2
u
r
∂
r
2
+
1
r
∂
u
r
∂
r
)
Φ
(
r
,
t
)
h
h
h
h
h
h
h
∂
2
Φ
(
r
,
t
)
∂
r
2
+
u
r
(
∂
u
t
∂
r
+
∂
u
t
∂
t
)
Φ
(
r
,
t
)
]
.
Before solving the above equation, one needs to relate the additional function in the modified equation to physical situation that takes place in the aquifers. Some deterministic models treat the properties of porous media as lumped parameters (essentially, as a black box), but this prevents the representation of heterogeneous hydraulic properties in the model. Heterogeneity or variability in aquifer properties is characteristic of all geologic systems and is now recognised as playing a key role in influencing groundwater flow and solute transport. Thus, it is often preferable to apply distributed-parameter models, which allow the representation of more realistic distributions of system properties.The lithology of most geological formations tends to vary significantly, both horizontally and vertically. Consequently, geological formations are seldom homogeneous. Figure1 is an example of layered heterogeneity.Figure 1
Example of heterogeneous karst aquifer illustrating the duality of recharge (allogenic versus autogenic), infiltration (point versus diffuse), and porosity/flow (conduits versus matrix) [12].Heterogeneity occurs not only in the way shown in the Figure1, however, individual layers may pinch out; their grain size may vary in horizontal direction, they may contain lenses of other grain sizes, or they may be discontinuous by faulting or scour-and-fill structures.The distribution of sedimentary facies controls the heterogeneity of hydrogeological properties of porous sedimentary aquifers at different scales. The arrangement of individual facies and their porosity and permeability determine the path of groundwater flow across sedimentary bodies. Therefore the capability to forecast hydrogeological heterogeneity due to facies changes helps to improve solutions of flow and diffusion problems in this kind of aquifer. When real aquifers are studied, it is impossible to model groundwater flow at a scale such that we can take into account the effects of fine-scale sedimentary heterogeneity; in fact this would require a precise knowledge of the sedimentary bodies that cannot be obtained from sparse data at some wells and this would be prohibitive for the required computing power. Therefore the fine scale heterogeneity is usually “up-scaled” and the heterogeneous real medium is substituted at a larger scale with an equivalent often anisotropic medium, whose parameters allow the reproduction of the average flow of the real heterogeneous sedimentary structure. In this paper the functionF
(
u
r
,
u
t
,
Φ
0
(
r
,
t
)
) will be considered to account for the effect of heterogeneity and variability of the geological formation system in which the groundwater flows.


## 3. Solutions of the Modified Groundwater Flow Equation

Numerical methods yield approximate solutions to the governing equation through the discretisation of space and time. Within the discretised problem domain, the variable internal properties, boundaries, and stresses of the system are approximated. Deterministic, distributed-parameter, and numerical models can relax the rigid idealised conditions of analytical models or lumped-parameter models, and they can therefore be more realistic and flexible for simulating fields conditions. Our next concern in this paper is to provide solution of the above equation. To achieve this we will make use of two techniques including: the green function and the variational iteration method. We will start with the variational iteration method.


### 3.1. Variational Iteration Method

The values of the variational iteration method and its applications for a range of categories of differentials equations can be viewed in [13–15].Following the work recently done by Theis in 1935 [16], in which they proposed an analytical solution to the standard version of the groundwater flow equation, this solution can be approximated as
(19)
Φ
0
(
r
,
t
)
=
Q
4
π
T
{
e
-
(
r
2
S
4
T
t
)
ln
⁡
[
1
+
α
4
T
t
r
2
S
]
}
,
where Q is the constant discharge rate, T is the transmissivity of the aquifer and Φ
0
(
r
,
t
) the piezometric head. It follows that the right side of (13a) and (13b) is known. On the basis of the above equation and knowing the function u
(
x
,
t
), one can derive a solution of (13a) and (13b) where the unknown is the functionΦ
1
(
r
,
t
). To make things simple, we put h
(
r
,
t
)
=
F
(
u
r
,
u
t
,
Φ
0
(
r
,
t
)
). And (13a) and (13b) become
(20)
∂
Φ
1
(
r
,
t
)
∂
t
=
K
S
0
[
1
r
∂
Φ
1
(
r
,
t
)
∂
r
+
∂
2
Φ
1
(
r
,
t
)
∂
r
2
]
+
h
(
r
,
t
)
.
To solve (20) by means of variational iteration method, we put (20) in the form
(21)
K
S
0
[
(
Φ
1
(
r
,
t
)
)
2
r
+
1
r
(
Φ
1
(
r
,
t
)
)
r
]
-
h
(
r
,
t
)
-
(
Φ
1
(
r
,
t
)
)
t
=
0
.
The correction functional for (21) can be approximately expressed as follows for this matter as
(22)
Φ
1
n
+
1
(
r
,
t
)
=
Φ
1
n
(
r
,
t
)
+
∫
0
t
λ
(
τ
)
{
K
S
0
[
(
Φ
n
1
(
r
,
τ
)
)
2
r
1
r
(
Φ
1
n
(
r
,
τ
)
)
r
∂
m
Φ
n
1
(
r
,
τ
)
∂
τ
m
h
h
h
h
h
h
h
h
h
h
h
h
h
+
1
r
(
Φ
1
n
(
r
,
τ
)
)
r
]
h
h
h
h
h
h
h
h
h
h
-
h
(
r
,
τ
)
-
∂
m
Φ
n
1
(
r
,
τ
)
∂
τ
m
}
d
τ
,
where λ is a general Lagrange multiplier [17], which can be recognized optimally by means of variation assumption [17–19], here (
Φ
1
(
r
,
t
)
)
ˇ
2
r, (
Φ
1
(
r
,
t
)
)
ˇ
r, and h
(
r
,
τ
)
ˇ are considered as constrained variations. Making the above functional stationary
(23)
δ
Φ
1
,
n
+
1
(
r
,
t
)
=
δ
Φ
1
n
(
r
,
t
)
+
δ
∫
0
t
λ
(
τ
)
{
∂
m
Φ
n
1
(
r
,
τ
)
∂
τ
m
}
d
τ
.
Capitulates the next Lagrange multipliers, giving up to the following Lagrange multipliers λ
=
-
1 for the case where m
=
1 and λ
=
t
-
τ for m
=
2. For these matter if m
=
1, we obtained the following iteration formula:
(24)
Φ
1
,
n
+
1
(
r
,
t
)
=
Φ
1
,
n
(
r
,
t
)
-
∫
0
t
{
K
S
0
[
(
Φ
1
,
n
(
r
,
τ
)
)
2
r
+
1
r
(
Φ
1
,
n
(
r
,
τ
)
)
r
]
h
h
h
h
h
-
h
(
r
,
τ
)
-
(
Φ
1
,
n
(
r
,
τ
)
)
τ
1
r
}
d
τ
.
Hence we commerce with
(25)
Φ
1,0
(
r
,
t
)
=
Φ
1
(
r
,
0
)
=
0
.
Means that before the water is pumped out from the borehole, the water level in the aquifer is the same and is considered here to be zero level.It is worth noting that if the zeroth componentΦ
0
(
r
,
t
) is defined, then the remaining components n≥ 1 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 Φ
(
r
,
t
) is approximated by the truncated series
(26)
Φ
1
N
(
r
,
t
)
=
∑
n
=
0
N
-
1
Φ
1
n
(
r
,
t
)
,
lim
⁡
N
→
∞
⁡
Φ
1
N
(
r
,
t
)
=
Φ
(
r
,
t
)
.
We follow next with the second component
(27)
Φ
1,1
(
r
,
t
)
=
∫
0
t
h
(
r
,
τ
)
d
τ
.
To calculate Φ
1,1
(
r
,
t
) we first need to define explicitly the function u
(
r
,
t
). The following function we define here does not actually have a physical meaning, but we use it as example.To make thing simple, we suppose thatu
t
=
1 and u
r
=
0.5 and the function h
(
r
,
t
) becomes
(28)
h
(
r
,
t
)
=
exp
⁡
[
-
r
2
S
4
t
T
]
(
Q
ln
⁡
[
1
+
4
t
T
α
r
2
S
]
-
1
)
×
(
-
16
Q
t
2
T
α
2
π
r
6
S
2
(
1
+
4
t
T
α
/
r
2
S
)
2
h
h
h
+
6
Q
t
α
π
r
2
T
(
1
+
4
t
T
α
/
r
2
S
)
h
h
h
+
2
Q
α
(
1
+
4
t
T
α
/
r
2
S
)
+
Q
r
2
S
2
ln
⁡
[
1
+
4
t
T
α
/
r
2
S
]
16
π
t
2
T
3
h
h
h
-
Q
S
ln
⁡
[
1
+
4
t
T
α
/
r
2
S
]
8
π
t
T
2
)
.
In this matter two components of the decomposition series were obtained of which Φ
(
r
,
t
) was evaluated to have the following expansion:
(29)
Φ
1
(
r
,
t
)
=
Φ
10
(
r
,
t
)
+
Φ
11
(
r
,
t
)
+
⋯


### 3.2. Green Function Methods

To solve (13a) and (13b), we go on to construct a suitable green’s function for this case in point. Let (
R
,
τ
1) be the green’s function to be constructed, where R
=
|
r
-
r
0
| and τ
1
=
|
t
-
t
0
|. G is chosen so as to satisfy homogeneous boundary conditions corresponding to the boundary conditions. It is important to notice that the homogeneous solution of (13a) and (13b) is similar to the diffusion equation if one replaces Φ
1
(
r
,
t
) by ψ
(
r
,
t
); therefore, the green function involved here is the green’s function for the diffusion equation. Since the aquifer is said to be infinite, the green function for flow equation for infinite aquifer is given by [20]
(30)
G
(
R
,
τ
1
)
=
4
π
T
S
0
(
1
2
π
τ
1
)
2
exp
⁡
[
-
T
2
R
2
4
S
0
2
τ
1
]
k
(
τ
1
)
.
Here the function k
(
τ
1
) is to be determined by using the boundary condition. The above equation satisfies an important integral property which is valid for n
=
2.Consider(31)
∫
G
(
R
,
τ
1
)
d
S
=
4
π
S
0
2
T
2
,
τ
1
>
0
.
This equation is an expression of groundwater flow. At a time and at a position, the piezometer is introduced in the borehole that taps the aquifer. The water that is pumped out from the aquifer through the borehole is migrating through the porous media, but in such a way that the total amount of water in the aquifer is reduced as time goes on if there is no recharge. Since (14) still holds, we can observe that
(32)
G
(
R
,
τ
1
)
⟶
4
π
S
0
2
T
2
δ
(
R
)
,
τ
1
⟶
0
.
In addition, the green’s function used for this purpose is a solution to the following equation:
(33)
∂
G
(
R
,
τ
1
∣
R
0
,
τ
10
)
∂
t
=
K
S
0
[
1
r
∂
G
(
R
,
τ
1
∣
R
0
,
τ
10
)
∂
r
+
∂
2
G
(
R
,
τ
1
∣
R
0
,
τ
10
)
∂
r
2
h
h
h
h
h
=
-
4
π
δ
(
R
)
δ
(
τ
1
)
1
r
∂
G
(
R
,
τ
1
R
0
,
τ
10
)
∂
r
]
.
The general solution of (13a) and (13b) can then be given as function of the green function as
(34)
Φ
1
(
r
,
t
)
=
∬
4
π
T
S
0
(
1
2
π
τ
1
)
2
e
x
p
[
-
T
2
R
2
4
S
0
2
τ
1
]
k
(
τ
1
)
h
(
R
,
τ
1
)
d
R
d
τ
1
.
Here
(35)
h
(
r
,
t
)
=
exp
⁡
[
-
r
2
S
4
t
T
]
(
Q
ln
⁡
[
1
+
4
t
T
α
r
2
S
]
-
1
)
×
(
-
16
Q
t
2
T
α
2
π
r
6
S
2
(
1
+
4
t
T
α
/
r
2
S
)
2
h
h
h
+
6
Q
t
α
π
r
2
T
(
1
+
4
t
T
α
/
r
2
S
)
h
h
h
+
2
Q
α
(
1
+
4
t
T
α
/
r
2
S
)
+
Q
r
2
S
2
ln
⁡
[
1
+
4
t
T
α
/
r
2
S
]
16
π
t
2
T
3
h
h
h
-
Q
S
ln
⁡
[
1
+
4
t
T
α
/
r
2
S
]
8
π
t
T
2
)
.
Due to lake of experimental data for this situation, no graphical representation will be presented in this paper. One needs to model the function of uncertainties introduced in this paper and use it for computational simulation and the analytical solution of the modified groundwater flow equation can then be compared with experimental data. Or one can from the standard solution measure the aquifer’s parameters and use it to determine the values of the function u
(
x
,
t
) and this is not done in this paper.


## 3.1. Variational Iteration Method

The values of the variational iteration method and its applications for a range of categories of differentials equations can be viewed in [13–15].Following the work recently done by Theis in 1935 [16], in which they proposed an analytical solution to the standard version of the groundwater flow equation, this solution can be approximated as
(19)
Φ
0
(
r
,
t
)
=
Q
4
π
T
{
e
-
(
r
2
S
4
T
t
)
ln
⁡
[
1
+
α
4
T
t
r
2
S
]
}
,
where Q is the constant discharge rate, T is the transmissivity of the aquifer and Φ
0
(
r
,
t
) the piezometric head. It follows that the right side of (13a) and (13b) is known. On the basis of the above equation and knowing the function u
(
x
,
t
), one can derive a solution of (13a) and (13b) where the unknown is the functionΦ
1
(
r
,
t
). To make things simple, we put h
(
r
,
t
)
=
F
(
u
r
,
u
t
,
Φ
0
(
r
,
t
)
). And (13a) and (13b) become
(20)
∂
Φ
1
(
r
,
t
)
∂
t
=
K
S
0
[
1
r
∂
Φ
1
(
r
,
t
)
∂
r
+
∂
2
Φ
1
(
r
,
t
)
∂
r
2
]
+
h
(
r
,
t
)
.
To solve (20) by means of variational iteration method, we put (20) in the form
(21)
K
S
0
[
(
Φ
1
(
r
,
t
)
)
2
r
+
1
r
(
Φ
1
(
r
,
t
)
)
r
]
-
h
(
r
,
t
)
-
(
Φ
1
(
r
,
t
)
)
t
=
0
.
The correction functional for (21) can be approximately expressed as follows for this matter as
(22)
Φ
1
n
+
1
(
r
,
t
)
=
Φ
1
n
(
r
,
t
)
+
∫
0
t
λ
(
τ
)
{
K
S
0
[
(
Φ
n
1
(
r
,
τ
)
)
2
r
1
r
(
Φ
1
n
(
r
,
τ
)
)
r
∂
m
Φ
n
1
(
r
,
τ
)
∂
τ
m
h
h
h
h
h
h
h
h
h
h
h
h
h
+
1
r
(
Φ
1
n
(
r
,
τ
)
)
r
]
h
h
h
h
h
h
h
h
h
h
-
h
(
r
,
τ
)
-
∂
m
Φ
n
1
(
r
,
τ
)
∂
τ
m
}
d
τ
,
where λ is a general Lagrange multiplier [17], which can be recognized optimally by means of variation assumption [17–19], here (
Φ
1
(
r
,
t
)
)
ˇ
2
r, (
Φ
1
(
r
,
t
)
)
ˇ
r, and h
(
r
,
τ
)
ˇ are considered as constrained variations. Making the above functional stationary
(23)
δ
Φ
1
,
n
+
1
(
r
,
t
)
=
δ
Φ
1
n
(
r
,
t
)
+
δ
∫
0
t
λ
(
τ
)
{
∂
m
Φ
n
1
(
r
,
τ
)
∂
τ
m
}
d
τ
.
Capitulates the next Lagrange multipliers, giving up to the following Lagrange multipliers λ
=
-
1 for the case where m
=
1 and λ
=
t
-
τ for m
=
2. For these matter if m
=
1, we obtained the following iteration formula:
(24)
Φ
1
,
n
+
1
(
r
,
t
)
=
Φ
1
,
n
(
r
,
t
)
-
∫
0
t
{
K
S
0
[
(
Φ
1
,
n
(
r
,
τ
)
)
2
r
+
1
r
(
Φ
1
,
n
(
r
,
τ
)
)
r
]
h
h
h
h
h
-
h
(
r
,
τ
)
-
(
Φ
1
,
n
(
r
,
τ
)
)
τ
1
r
}
d
τ
.
Hence we commerce with
(25)
Φ
1,0
(
r
,
t
)
=
Φ
1
(
r
,
0
)
=
0
.
Means that before the water is pumped out from the borehole, the water level in the aquifer is the same and is considered here to be zero level.It is worth noting that if the zeroth componentΦ
0
(
r
,
t
) is defined, then the remaining components n≥ 1 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 Φ
(
r
,
t
) is approximated by the truncated series
(26)
Φ
1
N
(
r
,
t
)
=
∑
n
=
0
N
-
1
Φ
1
n
(
r
,
t
)
,
lim
⁡
N
→
∞
⁡
Φ
1
N
(
r
,
t
)
=
Φ
(
r
,
t
)
.
We follow next with the second component
(27)
Φ
1,1
(
r
,
t
)
=
∫
0
t
h
(
r
,
τ
)
d
τ
.
To calculate Φ
1,1
(
r
,
t
) we first need to define explicitly the function u
(
r
,
t
). The following function we define here does not actually have a physical meaning, but we use it as example.To make thing simple, we suppose thatu
t
=
1 and u
r
=
0.5 and the function h
(
r
,
t
) becomes
(28)
h
(
r
,
t
)
=
exp
⁡
[
-
r
2
S
4
t
T
]
(
Q
ln
⁡
[
1
+
4
t
T
α
r
2
S
]
-
1
)
×
(
-
16
Q
t
2
T
α
2
π
r
6
S
2
(
1
+
4
t
T
α
/
r
2
S
)
2
h
h
h
+
6
Q
t
α
π
r
2
T
(
1
+
4
t
T
α
/
r
2
S
)
h
h
h
+
2
Q
α
(
1
+
4
t
T
α
/
r
2
S
)
+
Q
r
2
S
2
ln
⁡
[
1
+
4
t
T
α
/
r
2
S
]
16
π
t
2
T
3
h
h
h
-
Q
S
ln
⁡
[
1
+
4
t
T
α
/
r
2
S
]
8
π
t
T
2
)
.
In this matter two components of the decomposition series were obtained of which Φ
(
r
,
t
) was evaluated to have the following expansion:
(29)
Φ
1
(
r
,
t
)
=
Φ
10
(
r
,
t
)
+
Φ
11
(
r
,
t
)
+
⋯


## 3.2. Green Function Methods

To solve (13a) and (13b), we go on to construct a suitable green’s function for this case in point. Let (
R
,
τ
1) be the green’s function to be constructed, where R
=
|
r
-
r
0
| and τ
1
=
|
t
-
t
0
|. G is chosen so as to satisfy homogeneous boundary conditions corresponding to the boundary conditions. It is important to notice that the homogeneous solution of (13a) and (13b) is similar to the diffusion equation if one replaces Φ
1
(
r
,
t
) by ψ
(
r
,
t
); therefore, the green function involved here is the green’s function for the diffusion equation. Since the aquifer is said to be infinite, the green function for flow equation for infinite aquifer is given by [20]
(30)
G
(
R
,
τ
1
)
=
4
π
T
S
0
(
1
2
π
τ
1
)
2
exp
⁡
[
-
T
2
R
2
4
S
0
2
τ
1
]
k
(
τ
1
)
.
Here the function k
(
τ
1
) is to be determined by using the boundary condition. The above equation satisfies an important integral property which is valid for n
=
2.Consider(31)
∫
G
(
R
,
τ
1
)
d
S
=
4
π
S
0
2
T
2
,
τ
1
>
0
.
This equation is an expression of groundwater flow. At a time and at a position, the piezometer is introduced in the borehole that taps the aquifer. The water that is pumped out from the aquifer through the borehole is migrating through the porous media, but in such a way that the total amount of water in the aquifer is reduced as time goes on if there is no recharge. Since (14) still holds, we can observe that
(32)
G
(
R
,
τ
1
)
⟶
4
π
S
0
2
T
2
δ
(
R
)
,
τ
1
⟶
0
.
In addition, the green’s function used for this purpose is a solution to the following equation:
(33)
∂
G
(
R
,
τ
1
∣
R
0
,
τ
10
)
∂
t
=
K
S
0
[
1
r
∂
G
(
R
,
τ
1
∣
R
0
,
τ
10
)
∂
r
+
∂
2
G
(
R
,
τ
1
∣
R
0
,
τ
10
)
∂
r
2
h
h
h
h
h
=
-
4
π
δ
(
R
)
δ
(
τ
1
)
1
r
∂
G
(
R
,
τ
1
R
0
,
τ
10
)
∂
r
]
.
The general solution of (13a) and (13b) can then be given as function of the green function as
(34)
Φ
1
(
r
,
t
)
=
∬
4
π
T
S
0
(
1
2
π
τ
1
)
2
e
x
p
[
-
T
2
R
2
4
S
0
2
τ
1
]
k
(
τ
1
)
h
(
R
,
τ
1
)
d
R
d
τ
1
.
Here
(35)
h
(
r
,
t
)
=
exp
⁡
[
-
r
2
S
4
t
T
]
(
Q
ln
⁡
[
1
+
4
t
T
α
r
2
S
]
-
1
)
×
(
-
16
Q
t
2
T
α
2
π
r
6
S
2
(
1
+
4
t
T
α
/
r
2
S
)
2
h
h
h
+
6
Q
t
α
π
r
2
T
(
1
+
4
t
T
α
/
r
2
S
)
h
h
h
+
2
Q
α
(
1
+
4
t
T
α
/
r
2
S
)
+
Q
r
2
S
2
ln
⁡
[
1
+
4
t
T
α
/
r
2
S
]
16
π
t
2
T
3
h
h
h
-
Q
S
ln
⁡
[
1
+
4
t
T
α
/
r
2
S
]
8
π
t
T
2
)
.
Due to lake of experimental data for this situation, no graphical representation will be presented in this paper. One needs to model the function of uncertainties introduced in this paper and use it for computational simulation and the analytical solution of the modified groundwater flow equation can then be compared with experimental data. Or one can from the standard solution measure the aquifer’s parameters and use it to determine the values of the function u
(
x
,
t
) and this is not done in this paper.


## 4. Conclusion

In this paper we modify the standard version of groundwater flow by replacing the standard derivative with Riemann-Liouville fractional derivatives approximations. The modified equations (13a) and (13b) differ from the standard form of groundwater flow equation in three properties. There is a new operator that takes into account the variation in piezometric head and uncertainties function; second, the “force” appears due to the coordinate dependence of uncertainty function; and finally, there is a derivative-free term that depends only on the uncertainties time function. The modified equation takes into account both the flow through the porous media and the effect of variability of the aquifer or the concept of heterogeneity of the aquifers [21]. The modified equation was solved via the green function technique and variational iteration method.



---

*Source: 498381-2014-05-08.xml*