INTRODUCTION
Groundwater is an important natural resource. Many agricultural, domestic and industrial water users rely on groundwater as the sole source of low cost high quality water. However, in recent years it has become that many human activities can have a negative impact of both the quantity and quality of groundwater resource. For examples, the depletion of groundwater resource by excessive pumping and the contamination of groundwater by waste disposal and other activities.
Therefore, increased interest has been shown in studying the response of aquifers
due to heavy pumping. One way to objectively assess the impact of existing or
proposed activities on groundwater quantity and quality is to design mathematical
model describing the groundwater flow. Many studies on ground water flow modeling
are available (Huyakorn et al., 1994; Poeter
and Hill, 1997; Pohll et al., 1999; Winter
and Tartakovsky, 2002; Cloot and Botha, 2008). In
this study the steps required for the construction and numerical implementation
of such a model are investigated. Ground water flow models are useful for studying
various types of flow behavior by examining different aquifer problems (Mercer
and Faust, 1980).
In developing a groundwater flow model a conceptual model of the physical processes
assumed of governing the behaviour of the groundwater system has to be constructed.
Due to the complexity of groundwater phenomena it is not possible to implement
a conceptual model derived from basic principles of physics in full and a number
of simplifying assumptions have to be made. These include the assumption that
the motion of the water through the soil obeys Darcy’s law. The set of
model equations forms the basis of most models of groundwater flow through a
saturated soil matrix. The formulation is taking into account for the effect
of deformation of the soil matrix as well the coupling between this deformation,
the pressure head and the moisture content of the soil. However these basic
equations are still neither analytically nor numerically tractable and further
simplification has to be made in order to yield qualitative results.
Generally speaking analytical solution can be obtained only in the cases where
very restrictive assumptions are made. For example the additional assumption
that the governing soil matrix does not deform gives rise to a single equation
governing the pressure evolution of the system and can be solved analytically
for the case of radial flow. The analytical solutions to this specific problem
were first suggested by Muskat (1973) and Theis
(1935) and serve as important guides for interpreting the results of numerical
solutions obtained for situations in which less restrictive assumptions are
made.
In order to construct a model of fluid flow coupled with aquifer deformation
use is made of Biot’s poroelasticity theory (Biot, 1941)
to obtain the stress tensor of the soil matrix. Within this framework the stress
tensor for the soil matrix is assumed to be homogeneous, isotropic and to obey
the linear Hook’s law. Experimental evidence however suggests that induce
large strains in the soil matrix, display plastic behaviour are subject to permanent
deformations.
The resulting equations describing the deformation induced effects due to pumping are then numerically implemented for a radially symmetric system in which an aquifer is assumed to be bounded above and below by aquitards. A finite element method is then used to obtain a numerical solution of the problem. A number of numerical experiments to determine the nature of the deformation and pressure profiles in the aquifer system during the initial stages are carried out.
Finally the possible influence of the boundary condition applied into the radial
deformation at the borehole is investigated and the results are compared with
earlier results obtained by Hsieh (1996).
BASIC SET OF EQUATION FOR DEFORMABLE GROUND WATER FLOW MODEL
A reduced system of equations modelling the water flow through a porous deformable
medium is given by Alam (2005):
where, e is the dilatational strain, p is the pressure of the fluid, ρ_{w} is the density of water, β is the isothermal compressibility of water, ε is the porosity of the medium, K is the hydraulic conductivity tensor, φ is the piezometric head of the water, σ the symmetric stress tensor, θ_{w} is the volumetric fraction of water, F_{w} the strength of any sources or sinks of the water with in Ω, v^{k} the velocity of the specie k relative to the coordinate system (x_{1}, x_{2}, x_{3}), S is the surface forces per unit area of the surface, B is the body forces per unit volume of the body,v is the velocity of the body.
Although, these equations are mathematically considerably simpler than the more exact equations, it would still be more difficult to solve them analytically. A computer program based on the PetrovGalerkin finite element approximation of the equations was therefore developed to solve the equations numerically.
NUMERICAL METHODS
The model equations were solved numerically using a code referring to 3 different techniques, namely: A PetrovGallerkin finite elements method is used for integrating the system in space, a finite difference scheme is used to model the advanced of the solution in time and those two steps are imbedded in an iterative process which linearises the system of discretised equations. Those different aspects of the numerical treatment of the system (1) will be briefly discussed in the next subsections.
The petrovgallerkin finite elements method: The finite elements methods
are well established and documented numerical methods to solve the ground water
flow equations (Pinder and Gray, 1977). In the case of
the simplified consolidation model (1), this method creates a set of two nonlinear
equations of the form:
where, the matrices,
and vector,
have the following expressions:
The time discretisation: In this work the time derivatives appearing
in system (2) were discretised according to the Newman method (Smith,
1982).
If Δt is the time increment and
the solution for field
at time t = t_{i+1} = t_{i}+Δt and t = t_{i}, respectively
we have the following approximation for the first derivatives:
or, otherwise written
where, λ is a weighting coefficient whose value is lying between 0 and 1.
(Note, however that, this numerical scheme is unconditionally stable in the range of parameter λ ε [1/2, 1], only).
In a similar way one can also get an approximate for the second derivatives:
Insertion of Eq. 4 and 5 into 2 yields
then an implicit numerical scheme for the displacements, u_{k}, k =
1,2,3 and the pressure, p
Theoretically, this system of equations, as it stands, is solvable. However, the nonlinearity present in the problem is a source of great numerical difficulties and therefore one performed additional operations aiming at linearising Eq. 6.
Linearising the numerical scheme: Nonlinearity arises in the system
of Eq. 6 into two places, namely: The components of the vector
are functions of the stresses which, in their turn, are linked to the displacements
through the highly nonlinear relation:
and the matrix
involves the porosity of the medium, θ_{w}, which, within our framework,
is a function of the displacements, u_{k}. The procedures applied to
bypass those nonlinearities are described hereafter.
Linearising the Components of :
In order to linearise the components of
the nonlinear Hooke curve θ_{w}(σ, p) is replaced by a piecewise
linear function made of a set of, say, N linear stretches on which a linear
Hooke law holds. So doing, a particular stretch, say k, is characterized by
its position
and its slope,.
Therefore, for a specific state of deformation (ε_{ij}, σ_{ij})
the linear relations (1), hold:
Combining relations Eq. 6 and 8 gives,
after some straightforward mathematical manipulations, explicit expressions
relating linearly the stresses, σ_{ij}, the pressure, p and the
strains, εR_{ij}, according to the form:
where,
are relations which identify with those of a classical linear Hooke law but for a material which would be material characterized by a shear modulus:
On the other hand,
account for the corrections due to the nonlinearity:
and are depending on the characteristics of the stretches under consideration, only.
Linearising the matrix components M_{ij}: The matrix
involves the porosity of the medium, θ_{w}, which, within our framework,
is a function of the displacements, u_{k}. Therefore the components
appearing on the left hand side of equation (4.8) reads:
and are, strictly speaking, quadratic expressions of the displacements.
To circumvent this problem the following iterative process is used:
Assume that the solution at time t^{n} is known. A first approximate
of the solution at time t^{n+1},,
is computed using the porosity relevant for the previous step, i.e.,.
In a second step, the value of
is updated by the relation .
And a new approximate solution for the pressure and the displacements is recalculated.
Further on, the procedure described in step two is repeated until convergence
of the fields
within a prescribed tolerance, δ:
AQUIFER OF HSIEH
The first step was taken in evaluating the numerical algorithms to develop
a computer program in Fortran and apply it to a hypothetical aquifer similar
to the one considered by Hsieh (1996). The aquifer is
100 m thick and confined by two less permeable layers. The upper layer is also
100 m thick, the lower one extends to a depths of 5000 m. The water table in
the aquifer coincides with the land surface, which is the top of the first layer.
A borehole that penetrates all three units fully is cased throughout both confining
layers but screened over the thickness of the aquifer. Groundwater is pumped
at a constant rate of 5x10^{2} from the borehole.
INITIAL AND BOUNDARY CONDITIONS
Since, the aquifer system is radially symmetric, the numerical model was restricted to a vertical crosssection of the aquifer that extended from the borehole, with a radius of 0.1 m to an outer boundary at a distance of 5 km and a depth of 5 km. The outer boundary was considered as impervious and subjected to a constant stresses. The bottom boundary also considered as impervious and rigid, i.e., the boundary was not subject to displacements. The reason for this choice of the two boundaries was to ensure that they are sufficiently distant from the borehole, so as not to influence the behaviour of the aquifer near the borehole. The boundary at the land surface was considered free of applied forces and deformable and subject to a constant piezometric head. No radial displacements were initially allowed at the borehole casing, while the vertical component of the boundary traction was assumed to remain constant. This allowed the aquifer matrix to move vertically along the screen. It was further assumed that the pumping caused a uniform flux of water across the entire thickness of the aquifer, but that no water entered the borehole across the casing in the two confining layers.
FINITE ELEMENT IMPLEMENTATION
The finite element mesh used for 5 km by 5 km domain consists of 30 columns and 80 rows of rectangular finite elements with variable sizes to prevent numerical oscillations, which required elements as thin as 0.01 m along the interfaces between the aquifer and the confining layers. Numerical simulation performed with this finite element mesh is done using a linear stressstrain relation.
RESULTS AND DISCUSSION
The general behaviour of the aquifer system with the linear stressstrain relation is observed. The effect of boundary conditions for the Hshieh aquifer is also examined.
Behaviour of the aquifer system with linear stressstrain relation:
Figure 1 shows that the pressure dropped rapidly close to
the borehole when pumping began. As shown in Fig. 2 and 3,
this caused the rock matrix to deform in both horizontal and vertical directions.

Fig 1: 
Evolution of the piezometric pressure across the centre plane
of the aquifer at a few times after pumping the borehole began 

Fig 2: 
Evolution of the radial displacement (u_{r}) across
the centre plane of the aquifer at a few times after pumping the borehole
began 

Fig 3: 
Evolution of the vertical displacement (u_{z}) across
the interfaces of the aquifer and two confining layers at a few times after
pumping the borehole began 
The deformations in the radial direction (u_{r}) at first increased
continuously near the borehole, until they reached a maximum value, not far
from the borehole, after which they began decreased continuously. Notice also,
how the position of the maximum deformation in the radial direction moved away
from the borehole with time. The vertical displacement, u_{z}, on the
other hand, were essentially restricted to the interfaces of the aquifer and
confining layers, where the elastic parameters are discontinuous.
From those numerical simulations interesting results followed: The first, is that the vertical displacements are essentially restricted to the interfaces of the aquifer and confining layers where the elastic parameters are discontinuous. The second is that the deformation causes the aquifer to contract towards the screened segment of the borehole, both horizontally and vertically, while the confining layers contract horizontally and extend vertically. This contraction of the aquifer rocks and the vertical extension of the rocks in the confining layers must obviously affect not only the porosity of the aquifer system, but also pressure distribution and thus the ability of the aquifer to yield water to the borehole. The fact that the extend of the deformation decreases with distance, suggests that the deformation of an aquifer during pumping, may be responsible for very puzzling observation that the storativity values, derived from conventional aquifer tests, tend to decrease with distance from the borehole that is pumped during such a test. This conclusion is supported by the fact that the simulated drawdowns in the deformed model are considerably lower than in a conventional model, especially near the borehole. There is thus more water available near the borehole in the deformed model, than in the conventional model. However, this excess of water disappears with distance from the borehole.
Effect of borehole boundary conditions: Up to know the borehole wall
was not allowed to deform as observed by Botha et al.
(1998) and therefore no radial displacements were allowed along the borehole
wall. However, this boundary condition may sound very restrictive and quite
artificial with respect to the field observations. Therefore, we investigate
the matter by considering an alternative boundary condition where radial displacements
can freely take place at the borehole namely:
Similar numerical simulations as those done for the rigid case:
where, these performed for situations involving linear deformations.
A comparison of the numerical results with respect to the two types of boundary conditions prescribed hereafter for the Hsieh aquifers.
Effect of boundary condition for the aquifer of hsieh: In the Hsieh
aquifer (Hsieh, 1996) the effect of the boundary condition
and
show a definite effect in the direct neighbourhood of this borehole, but this
difference in the deformation is quickly fading away with the distance and the
two curves corresponding to the two kinds of boundary condition are at sight
identical for distances larger than 5 m away from the borehole (Fig.
4). This observation is further sustained by Table 1,
which provides the numerical differences between the radial and vertical displacements
as well as the variation in pressure where these quantities realise their maximum
values.

Fig. 4: 
Effect of Boundary condition on (u_{r}) with respect
to the radial distance from the borehole (time 2 h) 

Fig. 5: 
Time evolution of u_{r} at the borehole wall using
the boundary condition
in the middle of the aquifer (z = 150 m) 
Table 1: 
Comparison between u_{r} for the boundary conditions 

Like shown in Fig. 5 describing the time evolution of the amplitude on u_{r} at the borehole for the boundary condition Du_{r}r = r_{b} = 0, the amplitude of u_{r} is quickly increasing with time when pumping starts and reaches a saturation after about 1 h the pumping operation started. A full description of the profile of u_{r} along the vertical direction of the borehole is given in Fig. 6.

Fig. 6: 
Evolution of u_{r} along the vertical axis z and
along the borehole (time 2 h) 

Fig. 7: 
Effect of Boundary condition u_{r} with respect to
the radial distance from the borehole (time 2 h) 
Table 2: 
Comparison between Δu_{r}, Δu_{z},
ΔP for the boundary conditions u_{r} = 0, Du_{r} =
0 

Effect of boundary condition for a fractured aquifer: Effect of boundary
conditions for the fracture aquifer was found similar to the effect depicted
by the Hsieh aquifer (Hsieh, 1996). In the fractured
aquifer the effect of boundary condition is also limited to a close neighbourhood
of the borehole and the difference in the deformations is quickly fading away
with the distance and the two curves corresponding to the two kinds of boundary
conditions are identical for distances larger than 5 m away from the borehole
(Fig. 7). Table 2 provides the numerical
differences between the radial, vertical displacements as well as the variation
in pressure where these quantities realise their maximum values.

Fig. 8: 
Time evolution of u_{r} at the borehole wall using
the boundary condition D_{r}u_{r}r_{h} = 0 in the
middle of the aquifer (z = 21 m) 

Fig. 9: 
Evolution of u_{r} along the vertical axis z and
along the borehole (time 2 h) 
A similar behaviour as the one observed for the time evolution of u_{r}
in the case of the Hsieh model (Hsieh, 1996) was also
obtained as shown in Fig. 8. A full description of the profile
of u_{r} along the vertical direction of the borehole is given in Fig.
9.
CONCLUSION
A ground water flow model has been developed. The pressure is examined and
it is found that the pressure dropped rapidly close to the borehole after pumping
began. Effect of two boundary conditions
and
for the Hsieh aquifer in the deformation is coincides quickly with distance
larger than 5 m away from the borehole. Using the same boundary conditions the
effect in the fractured aquifer was found similar to the effect for the Hsieh
aquifer. Numerical simulations using different boundary conditions at the borehole
were performed. The results of these computations showed clearly that the choice
made for the boundary condition at the pumping site just affect quantitatively
the drawdowns and deformation profiles up to 5 m away from the borehole wall.