FINITE VOLUME METHOD OF MODELLING TRANSIENT GROUNDWATER FLOW

In the field of computational fluid dynamics, the f inite volume method is dominant over other numerica l techniques like the finite difference and finite el ment methods because the underlying physical qu antities are conserved at the discrete level. In the presen t study, the finite volume method is used to solve an isotropic transient groundwater flow model to obtain hydrauli c heads and flow through an aquifer. The objective is to discuss the theory of finite volume method and its applications in groundwater flow modelling. To ach ieve this, an orthogonal grid with quadrilateral contro l v lumes has been used to simulate the model using m xed boundary conditions from Bwaise III, a Kampala Surb urb. Results show that flow occurs from regions of high hydraulic head to regions of low hydraulic head un til a steady head value is achieved.


INTRODUCTION
The groundwater resources of the earth have for a long time been subjected to degradation as a result of man's increasing utilization of natural resources and worldwide industrialization. Beginning in the 1960s, contaminated aquifers were cleaned up and protected from further degradation in various countries around the world because government agencies identified groundwater as a valuable and increasingly important water resource (Batu, 2005). During this time, it was found that mathematical groundwater flow and solute transport modelling could be used as an efficient and cost-effective tool in the investigation and management of groundwater resources. Since then, mathematical models of groundwater flow have been widely used for a variety of purposes ranging from water supply studies to designing contaminant cleanup. The availability of computers and the development of efficient computer programs to do the computations involved in the models have also led to an increase in the use of numerical mathematical models in the analysis of groundwater flow and contaminant transport problems.
Mathematical models are conceptual descriptions or approximations that describe the physical system using mathematical equations. They based on solving an equation (or a system of equations) that describe the physical phenomenon. Such equations are called governing equations of the specified phenomenon. For groundwater flow, the governing equations are Darcy's law and the principle of mass balance (conservation).
Darcy's law is an equation that describes the flow of a fluid through a porous medium. The law was formulated in 1856 by French engineer Henry Darcy while working on a project involving the use of sand to filter the water supply for the city of Dijon in France. From his experiments (Fitts, 2012;Freeze, 1994), Darcy observed that the rate of flow through a homogeneous sand column of constant cross-sectional area was proportional to both the cross-sectional area of the column and the defference in water level elevations at the inflow and outflow reservoirs of the column and inversely proportional to the length of the column. This equation is usually written as Equation (1)  In 3D, Darcy's law is given as Equation (2) (Fitts, 2012): where, K x , K y and K z are the hydraulic conductivity values in the x, y and z direction respectively. If the hydraulic conductivity is independent of the direction of measurement at a point in the porous medium, that is, K x = K y = K z , the medium is called isotropic at that point. If the hydraulic conductivity varies with the direction of measurement at a point in a porous medium, that is, K x ≠ K y ≠ K z , the medium is called anisotropic at that point. It should be noted that real geologic materials are never perfectly homogeneous (isotropic) but to ease calculations, it's often reasonable to assume that they are (Fitts, 2012).
The law of mass conservation or continuity principle, states that there can be no net change in the mass of a fluid contained in a small volume of an aquifer. Any change in mass flowing into the small volume of the aquifer must be balanced by a corresponding change in mass flux out of the volume, or a change in the mass stored in the volume, or both. The continuity equation is derived by considering a very small part of an aquifer called a control volume having the shape of a rectangular parallel-piped box of dimensions ∇x, ∇y, ∇z centered at some point P(x, y, z), as shown in Fig. 1.
The quantity of water in the control volume can change when groundwater enters or leaves the control volume through the sides. A mass balance is obtained on the water flowing in and out of this control volume.
That is Equation (3): and can be expressed as Equation (4) (Delleur, 2010): where, S s is the specific storage coefficient and q = (q x , q y , q z ).
Combining Darcy's law (2) and the continuity Equation 4 gives Equation (5): Which is most universal form of the saturated flow equation, allowing flow in all three directions, transient , heterogeneous conductivities (for example K x = f(x)) and anisotropic hydraulic conductivities. Other less general forms of the saturated flow equation can be derived from Equation (5) This can be simplified further by making the assumption that the porous medium is homogeneous and isotropic, that is, K x = K y = K z = K and Equation 6 becomes: If the flow is steady state , the right-hand side of Equation 5-7 all become zero, that is Equation (8-10): ( ) In many applications, groundwater is modelled as two-dimensional in the horizontal plane. This is because most aquifers have an aspect ratio like a thin pancake, with horizontal dimensions that are hundreds of times greater than their vertical thickness (Fitts, 2012). Also, the bulk of resistance encountered along a typical flow path is resistance to horizontal flow. Thus, the groundwater hydrologist can assume the aquifer to be of constant thickness b and the flow to be horizontal (in the x-y plane). Hence, the flow equation for an isotropic homogeneous porous medium is given as: where, S = bS s is the storativity. The product bK is called the transmissivity and Equation (11) is often written as Equation (12): The Groundwater flow models are used to calculate the rate and direction of movement of groundwater through aquifers and confining units in the subsurface. These calculations are referred to as simulations. The outputs from the model simulations are the hydraulic heads and flow rates which are in equilibrium with the hydrogeological conditions (aquifer boundaries, initial and transient conditions and sources or sinks) defined for the modelled area.

MATERIALS AND METHODS
In this study, we consider the second-order transient groundwater flow Equation (13): And carry out a finite-volume simulation of the problem based on the boundary and initial conditions from (Herzog, 2007) for Bwaise III parish in Kawempe Division, Kampala District, Uganda. From (Herzog, 2007), the study area is bordered to the north by Nabweru road, to the east by Bombo road, to the south by the Bwaise-Nsooba drainage channel and the west by the Nakamiro drainage channel, as shown in Fig. 2. It is 6.65 ha large.
The aquifer in the study area was divided into two layers: Top layer A and bottom layer B and the hydraulic conductivity varies remarkably in the study area because of area buildings. The conductivity of the first layer is shown in Fig. 3 and the value for layer B was set to 0:017 m/d (meters per day). The bottom of the aquifer was defined as impermeable. The aquifer has a depth of 15 m meeting the bedrock that is impermeable. A simplified cross-section illustrating the groundwater flow system is shown in Fig. 4. (Herzog, 2007) (Herzog, 2007) The main flow enters the system from the eastern and northern borders. The groundwater leaves the system in western and southern directions. This assumption is supported by the direction of flow of flood surface water after rainfall events. Since the northern and eastern sides of the study area experience inflow, the arcs delineating those boundaries were assigned to be specified head arcs. Specified head boundary conditions make it possible to adjust the head at the boundary. In the first attempt to create a simulation, the specified head was set at 5m below the ground surface, uniformly along the boundary. The head is assigned at the nodes and varies linearly over the connecting arc.

Fig. 2. Bwaise III study area
The drainage channels that form the western and southern boundaries will act as sinks, i.e., remove water from the aquifer, as long as the groundwater table is above the elevation of the drain. The drain will have no effect if the groundwater level falls below the bottom elevation of the drain. The rate of flow from the aquifer to the drain is proportional to the difference in height between the groundwater table and the drain bottom. The constant of this proportionality is the conductance of the fill material surrounding the drain. In the wet season, the runoff drainage channels tend to flow full, with the water levels exceeding the groundwater table. Thus in the wet months, the conductance of the drain is set very low, to simulate the absence of flow from the groundwater into the drain. This ensures that groundwater does not flow into the drain considering that flow from the drain into the groundwater is most likely (Herzog, 2007).
In this study, we assume an isotropic study area with hydraulic conductivity K given by the average value of the hydraulic conductivities of the two layers. Thus, the parameters for Bwaise III study area used for our simulation are: The Storativity S is given as Equation (15) and the transmissivity T as Equation (16) Thus, we carry out a finite volume simulation using quadrilateral control volumes and orthogonal mesh of the problem Equation (17): with boundary and initial conditions Equation (18): h(300, y, t) h(x,300, t) 12 n. h 0; x 0, y 0, h(x, y,0) 5 where, n is the outward unit normal to the boundary. This model describes transient flow in a twodimensional homogeneous isotropic conned aquifer of constant thickness b.

Finite Volume Discretization
In computational fluid dynamics, a numerical solution approach is used to solve the coupled and nonlinear set of equations that characterize fluid flow problems (Noorbehesht and Ghaseminejad, 2013). The dominant numerical technique in computational fluid dynamics is the finite volume method. The basic methodology of the method involves three steps: • The domain is subdivided into a number of finitesized sub domains called control volumes and each control volume is represented by a finite number of grid points (Causon et al., 2011) • Integration of the governing differential equation over each control volume and applying the divergence theorem • Consideration of a profile assumption for the dependent variable to approximate the derivative terms resulting in a set of algebraic equations, one for each control volume Theorem 3.1 (Divergence Theorem) Let V be a simply connected region in the xy-plane enclosed by a piecewise smooth curve ∂V. Let n be the unit outward-pointing normal to ∂V. Then: where, dV is the element of area and dX is the element of length.
Applying the methodology onto model Equation 17, we can average Equation (17) by integrating it over an arbitrary control volume V as Equation (19): Applying the divergence theorem yields: where, S is the surface of the control volume and n represents the outward unit normal to the surface. Using quadrilateral control volumes on a uniform cartesian mesh (Fig. 5), the surface integral in Equation (20) can be split into the sum of the four surface integrals over the cell faces S c (c = e, w, n, s) of the control volume, such that Equation 20 can be equivalently written in the form: The integrals appearing in Equation (21) can be approximated by the average values at the midpoints of the faces and at the center of the control volume Equation (22) (Schafer, 2006): Substituting in Equation (21) yields Equation (23): Using a first order forward difference in time and using the fact that S e = S w =∆y , S n = S s = ∆x and V = ∆x∆y, Equation (23) can be written as Equation (24) The main challenge of the FVM is the approximation of the gradients (or uxes) at the cell faces. The accuracy of a control volume discretization depends heavily on the approximation of the ux at the midpoint of the control volume faces and many methods have been proposed to approximate the gradient along a control volume surface for different computational fluid dynamics applications (Loudyi et al., 2007;Jayantha and Turner, 2001;2003). To calculate the gradients at the midpoint of the cell faces, an approximate distribution of properties between nodal points is used (Versteeg and Malalasekera, 1995). The simplest and most obvious technique is the Central Dierencing Scheme (CDS) which which assumes that h is a linear function between any two node points and a second order approximation for the gradients is given Equation (27) (Schafer, 2006): The time level at which these derivatives are computed determines whether the scheme is explicit (k), implicit (k+1) or Crank-Nicholson (mixture of both previous and new time levels). Using an implicit scheme and substituting (27) into (26), we get the following Equation (28) which can be simplified to: Using the (i, j) notation, Equation 29

JMSS
which is the FVM equation at the interior node P or 2 ≤ i, j≤ N-1. Note that this equation is common to both the cell-centered FDM and the FVM.
To complete the scheme (30) we needed to update the formula also for the boundary cell nodes i = 1, j = 1, i = N and j = N. These were derived by taking the boundary conditions (18) into account. We introduced ghost cells i = 0, j = 0, i = N + 1 and j = N + 1 which were located just outside the domain (Fig. 6). The boundary conditions were used to ll these cells with values h 0,j , h i,0 , h N+1,j and h i,N+1 , based on the values h i,j in the interior cells. The same sheme (30) was then used also for I, j = 1 and I, j = N.
Consider the Neumann boundary condition h n ∂ ∂ = 0 at x = 0 and y = 0. We formally extend the definition of the solution h for (x; y) <0, that is, outside the domain. At x = 0 Equation (31): Or h x (t,0,y) = 0. Now Equation (32): Dropping the O ((∆x) 3 ) term, we got an expression for h 0,j in terms of h 1,j Equation (33): Similarly h y (x, 0, t) = 0 ) h i,0 = h i,1 . The Dirichlet boundary conditions h (300, y, t) = h(x, 300, t) = 0:5 can be approximated to second order by taking the average of two cells to approximate the value in between Equation (34): O(( y) ) 2 Leading to the approximation Equation ( Similarly, h N+1,j = 24-h N,j . If we let ∆y = ∆x, then scheme (30) becomes Equation (36): Rearranging terms with k+1 on the left and terms with k on the right hand side gives Equation (37): The star in the middle of The solution vector h was given by Equation (39) and the right hand side vector b as Equation (40) where, k and k + 1 are time levels.

JMSS
The value of the parameters ∆x, ∆t, S and T were known as is the value of the hydraulic head, K i, j h . The set of equations were solved simultaneously at each time step, starting from a set of initial conditions where h i,j is known for all (i, j) and proceeding through the time steps, k = 1, 2, 3,…..

RESULTS
Presented here are the results from the simulation of the model. Figure 8-13 present results obtained when the western and southern boundaries are impermeable while Fig. 14-20 present the results when an outflow of 50m 3 /day is allowed through the western and southern boundaries.

DISCUSSION
The results presented show that flow occurs from regions of high hydraulic head to regions of low hydraulic head until a steady head value is achieved, as shown in Fig. 8, 12, 14 and 18. This agrees with Darcy's conclusion that hydraulic head decreases in the direction of flow. The would be infow northern and esatern borders would also act as outflow borders if they are at a lower hydraulic head than the other points in the aquifer. Thus the major determinant of groundwater flow direction is hydraulic gradient.
When the western and southern borders were impermeable, a steady head value equal to the Dirichlet boundary condition at the northern and eastern borders is achieved. The rate at which points in the aquifer attain the steady head value depends on how close they are to the inflow borders. Points nearer the inflow borders attain steday state faster than those farther away from the inflow borders. When an outflow of 50 m 3 /day was allowed through the western and southern borders, water flows in or out of the aquifer until a steady head value, which in this case is at some points below the Dirichlet head boundary value at the northern and eastern borders, is attained. The closer a point is to the outflow borders, the lower its water level is below the steady Dirichlet head value of 12 m.

CONCLUSION
An orthogonal grid finite volume scheme applied to an isotropic transient groundwater flow model has been described in this study. We have used quadrilateral control volumes with nodes at the centers of the control volume and assumed that h varies linearly between any two nodal points. We have also seen that in the FVM scheme, conservation is guaranteed for each cell and this ensures that both local and global conservation are guaranteed no matter how coarse the mesh. A fully Science Publications JMSS implicit scheme has been used to approximate the derivatives. The spatial truncation error is O((∆x) 2 ) and the temporal truncation error is O(∆t), that is, the fully implicit scheme is second order accurate in space and first-order accurate in time. The scheme is also unconditionally stable. We have used direct methods to solve the discretized system.
It has been observed that water flows from regions at higher hydraulic head to regions at lower hydraulic heads. More accurate solutions would be obtained when a much finer mesh is used. The results obtained from our simulations also seem to agree with the simulations obtained when we use the finite element PDEtool in MATLAB.
In the study, we assumed that the study area was isotropic which in reality is not the case. We also assumed a fixed h value at the inflow borders which in nature is almost impossible. The inflow into our study area depends on whether its rainy season or dry season and there are always variations in the inflow of water depending on the season of the year. Further research would focus on using an anisotropic model and also using time dependent boundary conditions. An error analysis of our model would also be interesting to look at in future work.

ACKNOWLEDGEMENT
The researcher express their gratitude for the nancial support from the International Science Program (ISP) based at Uppsala University in Sweden.