Modeling an Aquifer: Numerical Solution to the Groundwater Flow Equation

We present a model of groundwater dynamics under stationary flow and governed by Darcy's Law of water motion through porous media, we apply it to study a 2D aquifer with water table of constant slope comprised of an homogeneous and isotropic media, the more realistic case of an homogeneous anisotropic soil is also considered. Taking into account some geophysical parameters we develop a computational routine, in the Finite Difference Method, that solves the resulting elliptic partial equation, both in a homogeneous isotropic and homogeneous anisotropic media. After calibration of the numerical model, this routine is used to begin a study of the Ayamonte-Huelva aquifer in Spain, a modest analysis of the system is given, we compute the average discharge vector as well as its root mean square as a first predictive approximation of the flux in this system, providing us a signal of the location of best exploitation; long term goal is to develop a complete computational tool for the analysis of groundwater dynamics.


I. INTRODUCTION
One of the most powerful tools to advance theoretical and practical knowledge in the characterization of water flow through ground porous media are computational numerical models, which must always be compared and calibrated with experimental measurements [1,2]. In the practice it is common to carry out these kind of studies with a wide combination of geophysical methods [3]. Physical models are often given in terms of partial differential equations, which for the specific case of groundwater dynamics they turn out to be analytically unsolvable when one tries to find a solution over realistic domains and conditions. The numerical solution, in the bast majority of the occasions, is used to predict the adequate management of hydric resources as can be seen in [4], where they present a numerical study of the Puebla Valley Aquifer in Puebla, Mexico.
Although actual aquifers are inhomogeneous and anisotropic, in order to construct a model, they are usually decompose in a great number of elements since it is feasible to solve certain problems assuming that each element behaves in its neighborhood as a small ideal aquifer [5]. In past years, there has been a lot of activity towards enhancing the way we model ground heterogeneity using geostatistical data and digital models or by the implementation of genetic algorithms [6]. Additionaly, there has been intensive work in the line of finding analytical solutions via new definitions of the derivative operator with fractional order, like the Caputo-Fabrizio derivative, these kind of outlines take into account the aquifer heterogeneity characterizing it as a scale problem. See for example [7], where they modify the groundwater flow equation replacing the conventional time derivative appearing in the transient term by the fractional derivative, thus allowing the description of material diffusion at different scales. This paper is organized as follows, in section II we give a quick review of the formulation of the Finite Difference Method and its numerical implementation, we describe briefly the key ideas of three of the main iterative schemes of solution. In section III we present the groundwater flow model, based on Darcy's Law, which is a good approximation to the dynamics of fluids in porous media provided they fulfill certain conditions. In section IV we present the results of implementing the iterative method to this problem and a brief discussion is given. Finally, in section V we give final remarks and conclusions along with future work perspectives.

II. FINITE DIFFERENCE METHOD AND THE ITERATIVE FORMULATION
The Finite Diference Method essentially transforms a differential equation into a system of algebraic equations by means of a spatial discretization of a physical problem's domain. This is performed considering a finite set of points on a rectangle, called grid, as shown in figure 1. Let h(x, y) be a dependent variable for which we must solve the problem, we will use the central points scheme, so we fix our atention into the point labeled with indexes (i, j) i.e., h i,j . Each index labeled point is called a node, so that in the FDM the derivatives are approximated using the Taylor series and computed in terms of nodal points. In order to show how this approximation is done let us focus on a single variable function f , so according to figure 2 if we Taylor expand this function to second order at points x = a + ∆x and x = a − ∆x we get substracting this two expressions and solving for f , ignoring higher order terms, we get where it is obvious the reason for the name central points scheme (or central differences scheme) for the procedure, the quantities for which one requires to solve are given in terms of the extreme points at each side. With this, equation (1) allows us to approximate numerically the first partial derivatives as Similarly, one can sum the Taylor expansions above and solve for the second derivative at x = a, f (a), from which the numerical approximations of the second partial derivatives are Let us say for the moment that we are interested in solving the Laplace equation in 2-dimensions therefore we replace the second partial derivatives for their Finite Difference numerical equivalents, given by eqs. (4) and (5), setting ∆x = ∆y = ∆ for simplicity, we arrive at which is the finite difference formulation of the Laplace equation, coincidentally as will be explained later this very equation corresponds to that of a steady groundwater flow in an homogeneous and isotropic aquifer. Along with appropriate boundary values one can compute quantities in the inner part of a grid as exemplified in figure 3 in order to solve the problem, in this case we have a set of boundary values on all the external edges of the grid, thus for this case we get which is now an algebraic system of equations instead of the original partial differential equation, once we solve this system we have an approximate set of values that describe the behavior of the function h(x, y) over the solution domain.
In practice one make a very large partition of the domain, since it is advisable that the partition size ∆ should be in the interval (0, 1) in order to satisfy convergence criteria and to guarantee a good physical approximation, therefore one ends up with a very large system of equations which is very hard if not feasible to handle analitically, so we are in need to implement an iterative numerical scheme of solution. In order to do so we solve eq. (6) in terms of h i,j , and with it we compute h(x, y) iteratively all along the grid, thus we obtain "test" values from which we construct new ones in every iteration, improving the quality of them on each cycle and tending to get a solution in a prescribed and physically acceptable numerical error tolerance. The quantity in eq. (7) is often refer to as the five-point operator in the literature, this corresponds to assign at the point of interest the mean value of the four nearest points to it in the nodal array, that is its first neighbors mean value.
There are various iterative methods to solve a system of linear equations, we will refer here to three of them: Jacobi iteration, Gauss-Seidel iteration and Successive Over Relaxation, although we present results only for the second one. These methods in general are the most efficient to solve equations containing a great number of unknowns since they get rid of the need to store data [8].

A. Jacobi Iteration
We can simply describe the Jacobi procedure as follows, to each interior point in figure 3 we apply eq. (7), this method does not uses a specific order to compute the h i,j values, we use the expresion where m labels the previous values to the iteration that it is been computed and labeled by m + 1. First one needs to propose a set of initial values, generally by guess but physically based on the boundary values of the problem. Following the example on figure 3 one is led to propose as an initial solution from which after three iterations we get one continues with the iterations until some physically based prescribed error tolerance is met.

B. Gauss-Seidel Iteration
Unlike the method described above Gauss-Seidel iteration is performed in an orderly fashion, we compute across the grid from left to right and downwards line per line so that on each node we can used the newest values to compute the next one during the same iteration, in this way we accelerate the convergence of the computation. The iterative expression one uses is the movement with which we traverse the grid is totally similar to that of reading a book's page. The fact we use the new available values makes possible that this method has better performance than the simple Jacobi iteration [9], as will be discussed in section IV.

C. Successive Over Relaxation
To implement the SOR method we must define the residue, c, this quantity is given by the change between two successive iterations in the Gauss-Seidel method, thus in every successive Gauss-Seidel iteration the residue softly vanishes on each node with the use of the new computed values, therefore converging more and more rapidly. In this method the residue is multiply by a relaxation factor ω with ω ≥ 1 such that the newly computed value is given by h m+1 in practice this means that in the SOR formula we substitute the definition of the residue along with the Gauss-Seidel expression for h m+1 We say then that one works in the over relaxation scheme since we add more residue to the value at hand, if on the contrary we have 0 ≤ ω ≤ 1 the updated value of h m+1 i,j is said to be under relaxed and it corresponds to a weighted average of the values used to compute the new one. Under relaxation is generally employed in order to make a non-convergent system converge by damping the oscillations in the computed values [9].

III. GROUNDWATER FLOW MODEL
In order to obtain the mathematical model of groundwater dynamics, we assume that water, or any other fluid under study, obeys Darcy's law of flow through porous media, that is to say that the fluid is incompressible, that its Reynolds number N R < 1, experimentally it has been determined that any fluid with a Reynolds number less than 1 satisfies Darcy's law and also that it is a good approximation for fluids with 1 ≤ N R < 10, and additionally one requires the flow to be laminar. In its simple and modern form Darcy's law relates q the specific discharge vector [m/s], σ the hydraulic conductivity tensor [m/s] and h(r) a scalar field standing for the total head [m] by means of where the negative sign indicates that the flow of water is in the direction of decreasing head. Assuming that our coordinate axes coincide with the principal axes of the conductivity tensor then for an isotropic homogeneous medium σ has a diagonal form with three identical constant elements and Darcy's law simplifies to q = −σ∇h, nevertheless the more realistic configuration involves an inhomogeneous and anisotropic hydraulic conductivity, in this cases this tensor can be represented as a 3 × 3 diagonal matrix where each component has a local value depending on the position. The equation governing groundwater flow can be obtained by the Control Volume approach [10] and it states that ∇ · q = S s ∂h ∂t + W , where S s is the specific storage and W represents flow in or out of the control volume, these two terms are usually written in a single function f (r, t) called infiltration and represents the presence of sources or wells, using eq. (9) this is written as which represents groundwater flow in an inhomogeneous anisotropic aquifer, this equation is solved in order to obtain the potential head scalar field h, from which we determine the rates and directions of flux from the specific discharge vector. In the case of steady flow in a homogeneous isotropic and confined aquifer, the continuity or conservation law assures that eq. (10) reduces to the Laplace equation.

FIG. 4:
Sketch of the mathematical model and boundary conditions for the regional groundwater system studied by Toth.
In the following we restrict ourselves to the case of a 2-dimensional confined aquifer and to the geometry proposed by Toth [11,12], where he studies a confined aquifer formed of a homogeneous, isotropic, porous medium with an underlaying layer of impermeable rock and a water table of linear slope. The system consists of a small watershed bounded on one side by a topographic high and on the other side by a mayor stream, both boundaries mark regional groundwater divides and they are treated mathematically as no-flow boundaries as well as the lower boundary, the upper boundary of the mathematical model is the horizontal line at height y = y 0 although the water table of the system lies above this line and the head along this boundary is assume to variate linearly as h(x, y 0 ) = cx + y 0 , where c is the slope of the water table. Thus, the domain of the problem is an approximation to the actual shape of a saturated flow region, in figure 4 we depict the cross section of the watershed along with the boundary conditions. The solution for this system in the simplified setup described above is given by cos (2m + 1) πx s cosh (2m + 1) πy needles to say that the homogeneity and isotropy simplifications are not compulsory to obtain a numerical solution.
In figure 5 we depict the profile of Toth's analytical solution.
It is in order to say that the study of analytical solutions for steady and transient flow, though maybe non totally realistic, can give us a great deal of insight of an aquifer and groundwater flow systems. For instance the stationary solution can be used to determine the response time of a system under transient flow, one of the methods to compute this time consists in solving both equations, for steady and transient flow, with which the response time will be the amount of time required for the transient solution to approach to the stationary one within a predefined tolerance [13][14][15]. Measure of this time scale is of critical importance since it can help us to define whether or not it is appropriate to use a transient flow model, if the time scale of interest is greater than the time response of the system then it is adequate to use a steady-state mathematical model and inversely [16]. Field work does not allow much room for error when establishing the computational model, in most cases is too expensive to redo the model should second guesses or doubts appear, at this point some simple analytical solutions, like those for one dimensional or radial flow, can provide us with a clear idea of the fundamental behavior of more complex or real world groundwater flow systems, thus serving as a tool both to decide between a steady or a transient model and to fine tune the parameters [16].
Moreover problems susceptible of being analytically solved are still the battleground to probe new predictive techniques, as the mean action integral recently proposed in [17]. Also, in the practice we are presented with the problem of sustainable exploitation, thus one is interested in the stationary solution and the response time of the system to the end of predicting the time it would take an aquifer to reach equilibrium once pumping has been introduced, specially if we want the pumping to continue essentially in an indefinite form, this time window in much cases ranges the millennia scale [18].

IV. NUMERICAL IMPLEMENTATION
Making use of the results of section II we will rewrite the general flow equation for a 2-dimensional aquifer as a purely algebraic discrete expression as follows, let us consider first the general case of an inhomogeneous anisotropic confined aquifer in steady flow, thus the right hand side of eq. (10) vanishes, also for the moment we leave the components of the hydraulic conductivity tensor unfixed, upon using eqs. (2)-(5) and solving for h i,j we get where we have set ∆x = ∆y = ∆ for simplicity.
To compute the solution of this system we use the Gauss-Seidel iteration method since it will converge faster than the simple Jacobi iteration, also given the typical values of h and the components of the σ tensor the algebraic equations represented by (12) are well posed and behaved, making unnecessary to use the Under Relaxation method, on the other hand, also in terms of the rate of convergence, according to these facts, one does not gain much if the Over Relaxation method were to be used.
We implement the numerical solution using Python 3.6, recalling we use a second order approximation let us say that we take ∆ = 0.1, therefore our approximation gives us a numerical result accurate to at most O(2) ∼ 10 −2 which means according to [19] that our calculations have three significant figures of accuracy, therefore we are forced to implement the iterative scheme within a loop that gradually rises the level of accuracy of the results until a prespecified error tolerance, e s , is achieved and with this the desired number of significant figures is reached.
For the aquifer setup described at the end of the previous section we have σ = constant, thus eq. (12) reduces to which is nothing else than the mean value of the four neighboring points of h i,j . Solving these numerically, we get the profile shown in figure 6, where we depict the head scalar field we obtain with our code, also figure 7 shows the super position of the analytical and numerical profiles, as we can see they are in great agreement showing overlapping in most of the domain region except in the far edges, where this behavior is to be expected given the numerical approximation scheme and its unavoidable natural error. In figure 8 we show the convergence curve of both Jacobi and Gauss-Seidel numerical solution methods for this homogeneous isotropic aquifer system, where we can see that Gauss-Seidel method reaches the prespecified error tolerance e s at approximately half the number of iterations with respect to the Jacobi method. This is no surprise if we recall the definition of the asymptotic rate of convergence of an iterative procedure, given by where B is the iteration matrix of the procedure and ρ is its spectral radius, then since ρ (B GS ) = ρ (B J ) 2 for these two methods [20], we have as a general result that λ GS /λ J = 2. This implies that in general Gauss-Seidel procedure converges twice as fast as Jacobi procedure, therefore in order to reach certain error tolerance it requires half the iterations with respect to the later.
These results make us confident to apply our computational code to a 2D-cross section of an homogeneous anisotropic aquifer within the same approximations as the problem described above. We take for instance the Ayamonte-Huelva aquifer, localized in the Spanish province of Huelva in Andalucía County, whose hydrogeological characteristics suggest the presence of materials such as sand, sandstone, clay, slate, gravel, lime and chalk [21,22]. As stated before the hydraulic conductivities of the soil and rocks play a main role in the dynamics of the system since conductivity depends on physical factors as the porosity, size and form of the grain, as well as the geometrical distribution of the particles in the media and the physical properties of the fluid, generally clay materials exhibit low values of hydraulic conductivity while sandstone and gravel have high values. The conductivity of saturated regions can be determined by various techniques including mathematical calculations, laboratory methods, tracer tests, well tests, auger hole tests and pumping tests of wells. Here we will used the representative values of vertical and horizontal hydraulic conductivity, σ y and σ x in our convention, for the rock types present in the Ayamonte-Huelva aquifer [23], they are 0.003721755 m/s and 0.003721761 m/s respectively, we computed these values in a rough approximation as the weighted average of the conductivity of each material composing the soil in the aquifer mentioned before. Under these homogeneity and anisotropy conditions eq. (12) is written as which in resemblance to the expression for the isotropic case this last equation amounts to a weighted mean of the first neighbors of h i,j , in figures 9 and 10 we depict the profiles of the numerical solution for this case, within the same boundary conditions as before, as well as its superposition with the isotropic case fig. 6. As we can see from figures 9 and 10 heads for both cases, isotropic and anisotropic aquifer, are basically equivalent which is to be expected since the numerical values of the horizontal and vertical conductivities we computed are almost equal, also their significant figures are the same for both values in the numerical model due to the error tolerance we prescribed; additionally the mean deviation of both calculations is order 10 −9 , thus lies bellow the accuracy we work with. This tells us that in practice we can treat this aquifer as an isotropic homogeneous one within our approximation. In figure 11 we also show the convergence curves upon implementation of the Gauss-Seidel and Jacobi iterative numerical procedures. From fig. 6 we can appreciate the qualitative behavior of the flow, recalling from Darcy's law that the negative sign in eq. (9) means that the groundwater flows in the direction of head loss, also the exact stream lines and the specific discharge vector can be seen from fig. 12. It is easy to see that at the boundary x = 0 km water flows to the left from the water table and in an analog way along the line x = 10 km water goes right from the water table, towards the side boundaries of the aquifer. Along the boundary y = 0 km flow is basically downwards with a slightly sideway deviation to the vertical boundaries of the aquifer. The y = −10 km receives a great amount of water, since as can be seen from the directional field approximately a third of the stream lines ends in this boundary. Using the numerical results we compute the average discharge vector, which in the units we used represents also the velocity of the flux, resulting in q = [0.001406x − 0.010008ŷ] m/s, this tells us, with in this approximation, that the most adequate region the exploitation of this aquifer is along the boundary on the righthand side, also upon computation of the root mean square of q we know that in general water flows at v = 0.019362 m/s, this figure is important since it helps us to estimate in an approximate fashion the volume of water we can extract per unit time. We should state that while groundwater flow in aquifers can be modeled via specialized softwares, based on finite differences and directly oriented to the problem, as Visual MODFLOW and MODFLOW-SURFACT in our opinion they present certain limitations, for instance MODFLOW requires additional subroutines, not included as part of the suit, in order to discretize the data one feeds to the software (via text archives) in the form it needs; on the other hand Visual MODFLOW's major disadvantage resides in the numerical formulation, it is impossible to "interpolate" between cell values when there is no measured or fed value available, this software drys out the cell assigning to it a zero head value, therefore producing a computational error; another problem is that its numerical resolution allows only to work with regular aquifers, generally rectangular or prismatic, therefore to model irregular, anisotropic or heterogeneous aquifers via this software is very difficult. All this situations can be overcome with careful, though laborious, coding and detailed calibration of the computer model, as well as a sufficiently accurate and precise finite difference scheme [20].
Also we must emphasize the role computer groundwater models have gained in a wide variety of applications from critic environmental problems to direct industrial applications, consider for instance a contaminated aquifer by heavy metals for which having a 3-dimensional transport model in the case of heterogeneous flow is of vital importance in order to implement soil remediation strategies, these models are generally described in terms of the water balance equation [24], a totally analogous treatment to the flow equation derived from Darcy's law. In the same sense an adequate groundwater model is of the uppermost relevance in productive activity like the extraction of rocks from a quarry. Understanding of the interaction between a quarry and groundwater is fundamental since rock operators in general excavate as deeper as the deposit allows, given this situation the water table of the formation is usually reached and pumping systems are implemented in order to dewater the aquifer, thus depressing the phreatic level beneath the bottom of the pit. Currently there are works such as [25] where interaction between quarry and groundwater is predicted by means of mathematical indexes developed from a multifold perspective, the key idea of using an index is that it quantifies in a single number the information of a 3-dimensional computational model could provide us with.

V. CONCLUSIONS
We have studied the numerical solution of the steady groundwater flow equation obtained from application of Darcy's law for flux in porous media, the precise computation of the discrete version of such equation was given considering the general case of an inhomogeneous anisotropic confined aquifer with steady flow. We tested our second order approximation scheme by confronting numerical results with a known analytical solution. Application of this methodology in the case of the Ayamonte-Huelva aquifer showed that and its actual hydraulic parameters, within the second order finite difference approximation of the differential operators and with the tolerance we work with, this aquifer can be regarded as comprised of an homogenous and isotropic medium since consideration of anisotropy produces negligible changes in the numerical values of the hydraulic potential (with respect to the isotropic case) and thus in the generic behaviour of water flow. We should emphasized that this is valid for the vertical and horizontal hydraulic conductivity values, computed as the weighted mean of the conductivities reported in the literature for the components of the aquifer's soil. With this setup we were able to determine a first approximation of the optimal location for exploitation of the aquifer as well as an estimation of the flux. Further work is required, we are preparing a 3-dimensional fourth order code to model Ayamonte-Huelva aquifer which incorporates actual water table data from several locations, soil properties according to depth and including extraction wells and recharge sources also, thus requiring matching conditions depending on the depth; setting the way for the end purpose of developing and alternative computational toolkit for the analysis of groundwater flow, this will allow us to predict the behaviour of the system and to give policies for the sustainable exploitation and management of the aquifer.