2. Materials and methods
2.7 Groundwater flow model
2.7.1 General concepts of numerical modelling
Computer models are popular tools in groundwater management to evaluate the impact of any intervention on groundwater quantity and quality. In regard to the complex nature of subterraneous flow, the first approach is to develop a conclusive conceptual model of the hydrogeological processes within the model area. The second step is the transfer of this model into a mathematical model. The mathematical description is made by partially differential equations (PDE) and auxiliary conditions (ISTOK 1989) which are solved by numerical methods. With modern microcomputers it is possible to solve these equations in convenient time.
The governing partial differential equations used in mathematical models for groundwater flow is the Laplace’s equation for the saturated steady state condition (Eq.
2.8) and its derivative for the saturated transient state (Eq. 2.9):
h h h
K x K y K z 0
x x y y z z
⎛ ⎞
⎛ ⎞ ⎛ ⎞
⎜ ⎟
⎜ ⎟ ⎜ ⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎝ ⎠ ⎝ ⎠ ⎝ ⎠
∂ ∂ + ∂ ∂ + ∂ ∂ =
∂ ∂ ∂ ∂ ∂ ∂ (Eq. 2.8)
s th
h h h
K x K y K z S
x x y y z z
⎛ ⎞
⎛ ⎞ ⎛ ⎞
⎜ ⎟
⎜ ⎟ ⎜ ⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎝ ⎠ ⎝ ⎠ ⎝ ⎠
∂∂
∂ ∂ + ∂ ∂ + ∂ ∂ =
∂ ∂ ∂ ∂ ∂ ∂ (Eq. 2.9)
where, t = time;
h = hydraulic head;
Q = volumetric flux per unit volume representing source/sink terms;
Kx, Ky, Kz = hydraulic conductivity along the x,y,z axes which are assumed to be parallel to the major axes of hydraulic conductivity;
Ss = specific storage defined as the volume of water released from storage per unit change in head per unit volume of porous material.
Laplace’s equations combine Darcy’s law and the continuity equation into a single second-order partial differential equation (WANG and ANDERSON 1982). In mathematical correct terms it would be necessary to write the hydraulic conductivity K with negative signs to express that the flow is directed to the bottom. However, any conceptual mathematical description of a complex natural environment is necessarily a greatly simplified description of the reality (ISTOK 1989). The auxiliary conditions mentioned above are the initial conditions and the boundary conditions (MARSAL 1989). They will be highlighted in Chapter 8.
There exist several types of numerical models. The most popular among them are the:
• Finite difference models (FD models) and,
• Finite element models (FE models).
The first visual difference among these two methods is the way to discretise the model environment (see Fig. 2.9).
The FD method applies a regular element mesh with constant cell sizes. It works best for rectangular or prismatic aquifers of uniform composition (ISTOK 1989). Its advantage is the easily understandable mathematical formulations and it’s solving algorithms.
Especially the MODFLOW code developed by the USGS is popular and frequently applied (see e.g. CHIANG and KINZELBACH 2000; HARBAUGH et al. 2000).
Fig. 2.9: Some types of elements, based on the finite difference concept and the finite element concept (from ANDERSON and WOESSNER 1994).
The method of finite elements is principally used in mechanical engineering. Since the late 1970s applications for groundwater modelling were developed (WANG and ANDERSON 1982). Although the FE method is based on the same basic PDEs like FD it uses other mathematical ways to solve them.
In the FE method the defined domain (solution) of a PDE is replaced by a discrete matrix of points and each ablation by a differential expression which refers to the point matrix (MARSAL 1989). The domain can be structured in partial areas of almost arbitrary size.
These parts are called the finite elements. The shape of finite elements is mathematically flexible but most of the available software programs use the classic triangular shape (Fig. 2.10).
Fig. 2.10: The typical triangular element e. Each node is labelled (i, j, m) counter-clockwise and has its own x, y-coordinates marked by the specific footnote (modified from WANG and ANDERSON 1982).
The nodes of each element are very important. Mathematically, the element nodes are the homologues of the matrix points of the FD method (MARSAL 1989).
The typical solution method for the FE method is Galerkin’s method (e.g. WANG and ANDERSON 1982; ISTOK 1989).Galerkin uses the weighted residual principle to govern the PDEs for groundwater flow. It is an integral formulation which leads to a set of algebraic equations that can be solved for values in the field variable (piezometric head h and/or the concentration C for transport problems) at each node in the mesh.
First an approximate solution to the boundary or initial value problem is defined. When this approximate solution is introduced into the governing differential equation, an error the so-called residual occurs at each point in the problem domain. Then the weighted average of the residuals for each node in the finite element mesh is forced to equal zero.
Despite more complex mathematics and more computational demand the FE method offers some advantages especially for more complex models (ISTOK 1989):
1. Irregular or curved aquifer boundaries, anisotropic and heterogeneous aquifer properties and sloping soil and rock layers can easily be incorporated into the numerical model by refining the mesh and individual positioning of the elements.
2. The accuracy of solutions to groundwater flow and transport problems is very good.
3. Solutions of transport models are generally more accurate than solutions obtained by the FD method.
4. The numerical stability of equation solvers is generally better than in FD models.
2.7.2 Modelling standards of FEFLOW®
The regional groundwater modelling was realised with the modelling software package FEFLOW® 5.1 from WASY GmbH. Details about actual versions and computing capabilities can be found on the WASY homepage5. DIERSCH (2005) gives detailed information about all numerical methods, solvers and mathematical methods available in this software. FEFLOW® was chosen for different reasons:
• Due to the great size of the HVO some areas stay unknown while at other locations detailed research was made. The flexible mesh refinement due to the FE method is an advantage in this case. New field investigations (e.g. pumping tests) may deliver new information about some areas. The FE mesh can easily be modified and thus may develop together with the ongoing field research. However, practical modelling needs a number of reruns to find out optimal design variables (DIERSCH 2005).
• All spatial data can be produced externally by GIS software. The FEFLOW® interface accepts a great variety of data formats as input and also as output. Exchange of data is therefore easy and fast-forward.
• Besides the many features of FEFLOW® there is a wide choice between different equation solvers and adaptive alternative numerical solutions. Accelerated solvers allow the solution of large problems with either great number of nodes or complex structures (DIERSCH 2005).
All input data can be handled on element level. Transient and intransient data can be entered. Point data can be assigned to the mesh by geostatistical interpolation techniques such as Akima, kriging or inverse distance or simpler by polygon joining
5 URL: http://www.wasy.de/english/index.html
techniques (DIERSCH 2005). For the purpose of this thesis the direct transfer of data from GIS-shape-files on the FE nodes, the so-called joining method, was preferred. This method was applied for as many parameters as for e.g. recharge, conductivity and storage coefficient.
Topographic data was preprocessed with the geostatistical software SURFER. A variogram and a regular data grid were produced and then imported into FEFLOW® by a nearest neighbour function to assure that nodal elevation values are as close as possible to their real distribution.
The regional groundwater recharge is calculated externally. FEFLOW® does not contain a special recharge or evapotranspiration package as it is used for example in MODFLOW applications (e.g. HARBOUGH et al. 2000).
The regional model is realised in 3D with three layers. For the model saturated conditions with an unconfined and movable upper surface were selcted. The automatic time-step control by the 2nd order predictor-corrector scheme “forward Adam-Bashforth”
was chosen with an initial time step length of 0.001 d. Time step control was given without any further constraints. The Galerkin-based finite element method for unstructured meshes and the iterative PCG (preconditioned conjugate gradient) solver are chosen. Other convergence and error criteria are set to default.