## Numerical groundwater models 1551 Assumptions and boundary conditions

Complex and transient (non-steady) flow conditions require an approximate numerical solution to the groundwater flow equations. There are a variety of numerical solution methods available, including finite difference, finite element, finite volume and boundary element methods. All are general methods for solving partial differential equation, such as that based on Darcy's law used to describe groundwater flows. Similar solvers are now available in many general mathematics packages such as Matlab, Mathcad and Mathematica. The finite-difference method is the most straightforward and easily understood numerical method and is the basis for groundwater packages such as the widely used MODFLOW package, originally developed by the US Geological Survey but now also available with graphical interfaces from a variety of different vendors (such as Visual MODFLOW). MODFLOW, and other groundwater packages can be used to solve problems in one, two and three dimensions with a variety of boundary conditions and for both steady-state and transient flow conditions.

In the finite difference method, the groundwater flow domain is subdivided into blocks in a similar way to that used to derive the flow equations in Fig. 15.2 (e.g. Fig. 15.8). In solving the equations, however, it is more convenient to use a grid of nodes at the centre of each block. The nodes of the grid are the points at which the values of head are required or can be specified as boundary conditions for the flow domain. At each node, the h value is assumed to be representative of the block centred on that node. If the value of h at A is unknown, it can be calculated from considerations of the flow pattern around the block.

 C A i D a ; qda qeaj x - axâ€”- F e G x

Fig. 15.8 Definition diagram for a block-centred finite difference scheme. The block around the point A is indicated by dashed lines. Flows Q are to be calculated given heads at the centres of the neighbouring blocks, B, C, D and E.

Fig. 15.8 Definition diagram for a block-centred finite difference scheme. The block around the point A is indicated by dashed lines. Flows Q are to be calculated given heads at the centres of the neighbouring blocks, B, C, D and E.

Then the finite difference estimate of the flux from A to B is

To obtain a complete finite difference form for the full transient flow equation (15.9), we need to define finite difference approximations for differential terms of the form d / 9h\ 9x\ * bxj

centred around the point A. This can be achieved by using a Taylor series expansion of h around the point A. By including higher order terms in the expansion, we can get a closer approximation to the values of the differential terms. Most finite difference models use a second-order approximation for the spatial derivatives of the form (assuming that Kx is constant everywhere):

Similar expressions can be used for the differentials in the y and z directions. If Kx is not isotropic, then the change in the conductivities in the different directions can be included in formulating the flux terms (this is clearly easiest if the nodes are on the boundaries between regions of the flow domain of different conductivity). If the conductivities are not homogeneous because of changes in aquifer properties, lithology, occurrence of faults or other features, then different conductivities can be used in different elements of the discretisation of the flow domain.

For steady-state flows the spatial difference approximations of (15.20) will be sufficient to formulate a finite difference solution for all the nodes in the grid. Effectively, since the solution at A will depend on the solution at the surrounding points B, C, D, E, which are also initially unknown except at fixed head boundary points, the result is a set of linear simultaneous equations, including equations for the boundary conditions. Solution techniques for such cases are well known, and involve iterating the solution from an initial estimate until a final solution is reached.

For transient problems, we also need to define a finite difference estimate for the time differential 9h/dt in (15.9).

We also need to now index the nodal head values by a time step index, t = 1, 2, 3 T

for time steps of length At. There are two approaches to solving such time-stepping problems. The easiest to implement is an explicit solution. This makes the solution at time t a function only of the known values at the previous time step at t - At (Fig. 15.9a). In this case, given a pattern of heads for the first time step (the initial conditions) at time t = 0, the spatial differential terms are evaluated at time t = 0, and the solution at time t = At is obtained using the forward difference approximation to

Fig. 15.9 Schematic diagram for: (a) explicit; and (b) implicit time-stepping schemes for finite difference representation of the time derivative. Filled circles represent values that are known at time step tj. Open circles represent values to be calculated. Arrows indicate dependencies in the calculation of hA t=At.

the time derivative d hA _ hA,t=At â€” hA,t=0 (1521)

Then having solved the set of equations for time t = At, the process can be repeated for successive time steps. The calculations are relatively simple because, at each time step, the only unknown in each nodal equation is now hA t, all the other terms use the values of h at time t â€” At, which are known. This formulation of the time stepping is only first order correct and can be subject to stability problems if there are strong gradients in time or space, or rapid changes in boundary conditions, unless rather short time steps are used.

An implicit solution can also be used (Fig. 15.9b). In this case, at each time step the time differential is treated in a similar way but the spatial differentials are treated as a weighted sum of the values at the times t and t + At. Thus, in an implicit formulation for a typical node at position A:

d K dJh) = VKx d x \ x d x) x hD,t+At 2hA,t+At + hB,t+At (Ax)2

where ^ is a time-weighting coefficient. The values of h at time t are known at the start of the time step, but those at time t + At are not. This leads to a system of linear simultaneous equations in the unknown values of h at time t + At. Efficient techniques for solving such equations are well developed, but involve a computational cost greater than for an explicit solution for each time step. The advantage of using an implicit solution, however, is that it is much less likely to suffer from stability problems so that much longer time steps can be used. The overall computational cost may then be less. The USGS groundwater modelling package, MODFLOW, for example, uses an implicit finite difference approximation of the continuous differential flow equations for both confined and unconfined aquifer cases.

For both explicit and implicit finite difference approximations, the starting point at time t = 0 does depend on being able to specify a set of initial conditions. Accuracy in reproducing the behaviour of the real system will depend on accuracy in the initial conditions. In some cases it will be possible to interpolate the pattern of initial conditions from measured heads from observation wells. In other cases, a steady-state solution might be used as the starting point for a transient solution, or a 'run-in' period might be used to obtain a consistent initial condition before simulating a period of interest.