Subscribe Now Subscribe Today
Research Article
 

Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer



Md. Johurul Alam, Mohammed Ashaque Meah and Md. Shah Noor
 
Facebook Twitter Digg Reddit Linkedin StumbleUpon E-mail
ABSTRACT

A conceptual model of the physical processes assumed of governing the behaviour of the ground water system has constructed to develop a ground water flow model. The resulting equations describing the deformation induced effects due to pumping are 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 used to solve the model equations numerically. 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. Effect of boundary conditions for the Hsieh aquifer in the deformation is examined. Using the same boundary conditions the effect in the fractured aquifer was also examined and then compared to the effect for the Hsieh aquifer. Numerical simulations using different boundary conditions at the borehole were performed. The study found that the developed flow model gives satisfactory results with that of Hsieh.

Services
Related Articles in ASCI
Search in Google Scholar
View Citation
Report Citation

 
  How to cite this article:

Md. Johurul Alam, Mohammed Ashaque Meah and Md. Shah Noor, 2011. Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer. Asian Journal of Mathematics & Statistics, 4: 33-44.

DOI: 10.3923/ajms.2011.33.44

URL: https://scialert.net/abstract/?doi=ajms.2011.33.44
 
Received: August 15, 2010; Accepted: September 25, 2010; Published: November 08, 2010



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

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
(1)

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, Fw the strength of any sources or sinks of the water with in Ω, vk the velocity of the specie k relative to the co-ordinate system (x1, x2, x3), 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 Petrov-Galerkin 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 Petrov-Gallerkin 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 sub-sections.

The petrov-gallerkin 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 non-linear equations of the form:

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
(2)

where, the matrices, Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer and vector, Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer have the following expressions:

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
(3)

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 Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer the solution for field Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer at time t = ti+1 = ti+Δt and t = ti, respectively we have the following approximation for the first derivatives:

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer

or, otherwise written

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
(4)

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:

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
(5)

Insertion of Eq. 4 and 5 into 2 yields then an implicit numerical scheme for the displacements, uk, k = 1,2,3 and the pressure, p

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
(6)

Theoretically, this system of equations, as it stands, is solvable. However, the non-linearity 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: Non-linearity arises in the system of Eq. 6 into two places, namely: The components of the vector Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer are functions of the stresses which, in their turn, are linked to the displacements through the highly non-linear relation:

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
(7)

and the matrix Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer involves the porosity of the medium, θw, which, within our framework, is a function of the displacements, uk. The procedures applied to by-pass those non-linearities are described hereafter.

Linearising the Components of Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer: In order to linearise the components of Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer the non-linear 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 Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer and its slope,Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer. Therefore, for a specific state of deformation (εij, σij) the linear relations (1), hold:

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
(8)

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, εRij, according to the form:

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
(9)

where,

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer

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:

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
(10)

On the other hand,

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer

account for the corrections due to the non-linearity:

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
(11)

and are depending on the characteristics of the stretches under consideration, only.

Linearising the matrix components Mij: The matrix Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer involves the porosity of the medium, θw, which, within our framework, is a function of the displacements, uk. Therefore the components

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer

appearing on the left hand side of equation (4.8) reads:

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer

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 tn is known. A first approximate of the solution at time tn+1,Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer, is computed using the porosity relevant for the previous step, i.e.,Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer.

In a second step, the value of Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer is updated by the relation Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer.

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 Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer within a prescribed tolerance, δ:

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer

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 cross-section 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 stress-strain relation.

RESULTS AND DISCUSSION

The general behaviour of the aquifer system with the linear stress-strain relation is observed. The effect of boundary conditions for the Hshieh aquifer is also examined.

Behaviour of the aquifer system with linear stress-strain 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.

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
Fig 1: Evolution of the piezometric pressure across the centre plane of the aquifer at a few times after pumping the borehole began

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
Fig 2: Evolution of the radial displacement (ur) across the centre plane of the aquifer at a few times after pumping the borehole began

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
Fig 3: Evolution of the vertical displacement (uz) 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 (ur) 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, uz, 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:

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer

Similar numerical simulations as those done for the rigid case:

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer

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 Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer and Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer 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.

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
Fig. 4: Effect of Boundary condition on (ur) with respect to the radial distance from the borehole (time 2 h)

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
Fig. 5: Time evolution of ur at the borehole wall using the boundary condition Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer in the middle of the aquifer (z = 150 m)

Table 1: Comparison between ur for the boundary conditions
Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer

Like shown in Fig. 5 describing the time evolution of the amplitude on ur at the borehole for the boundary condition Dur|r = rb = 0, the amplitude of ur 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 ur along the vertical direction of the borehole is given in Fig. 6.

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
Fig. 6: Evolution of ur along the vertical axis z and along the borehole (time 2 h)

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
Fig. 7: Effect of Boundary condition ur with respect to the radial distance from the borehole (time 2 h)

Table 2: Comparison between Δur, Δuz, ΔP for the boundary conditions ur = 0, Dur = 0
Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer

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.

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
Fig. 8: Time evolution of ur at the borehole wall using the boundary condition Drur|rh = 0 in the middle of the aquifer (z = 21 m)

Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer
Fig. 9: Evolution of ur along the vertical axis z and along the borehole (time 2 h)

A similar behaviour as the one observed for the time evolution of ur in the case of the Hsieh model (Hsieh, 1996) was also obtained as shown in Fig. 8. A full description of the profile of ur 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 Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer and Image for - Numerical Modeling of Ground Water Flow and the Effect of Boundary Conditions for the Hsieh Aquifer 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.

REFERENCES

1:  Alam, M.J., 2005. Modelling the effect of soil deformations on pumped aquifers, the linear case. M.Sc. Thesis, University of Free State, South Africa.

2:  Biot, M.A., 1941. General theory of three-dimensional consolidation. J. Applied Phys., 12: 155-164.
Direct Link  |  

3:  Botha, J.F., J.P. Verwey, I. van der Voort, J.J.P. Vivier, J. Buys, W.P. Colliston and J.C. Loock, 1998. Karoo aquifers. Their geology, geometry and physical behaviour. WRC Report No. 487/1/ 98. Water Research Commission, Pretoria.

4:  Cloot, A. and J.F. Botha, 2008. A generalised groundwater flow equation using the concept of non-integer order derivatives. Water SA., 32: 1-8.
Direct Link  |  

5:  Hsieh, P.A., 1996. Deformation-induced changes in hydraulic head during groundwater withdrawal. Ground Water, 34: 1082-1089.
Direct Link  |  

6:  Huyakorn, P.S., S. Panday and Y.S. Wu, 1994. A three-dimensional multiphase flow model for assesing NAPL contamination in porous and fractured media, 1. formulation. J. Contaminant Hydrol., 6: 109-130.
CrossRef  |  

7:  Mercer, J.W. and C.R. Faust, 1980. Ground-water modeling: An overview. Ground Water, 18: 108-115.
Direct Link  |  

8:  Muskat, M., 1973. The Flow of Homogeneous Fluids through Porous Media. McGraw-Hill Co. Inc., USA., pp: 55-120, 621-667

9:  Pinder, G.F. and W.G. Gray, 1977. Finite Element Simulation in Suface and Subsurface Hydrology. Academic Press Inc., New York

10:  Poeter, E.P. and M.C. Hill, 1997. Inverse models: A necessary next step in ground-water modeling. Ground Water, 35: 250-260.
CrossRef  |  

11:  Pohll, G., A.E. Hassan, J.B. Chapman, C. Papelis and R. Andricevic, 1999. Modeling ground water flow and radioactive transport in a fractured aquifer. Ground Water, 37: 770-784.
Direct Link  |  

12:  Smith, I.M., 1982. Programming the Finite Element Method. John Wiley and Sons, London

13:  Theis, C.V., 1935. The relation between the lowering of the Piezometric surface and the rate and duration of discharge of a well using ground-water storage. Trans. Am. Geophys. Union, 16: 519-524.

14:  Winter, C.L. and D.M. Tartakovsky, 2002. Groundwater flow in heterogeneous composite aquifers. Water Res. Res., 38: 1-11.
Direct Link  |  

©  2021 Science Alert. All Rights Reserved