Numerical Model Setup

Once a conceptual model has been developed (Section 4) and the most suitable modelling approach has been determined to be numerical methods, the numerical code must be selected (Section 5) and the conceptual model has to be implemented into a numerical model (see Figure 2-1). This step of the modelling process is often referred to as "model setup" but is also known as "model development" or "model construction".

In essence, the numerical model setup represents the process of converting a qualitative conceptual model into a numerical model, i.e. a complex set of mathematical equations that can be solved numerically. In order to solve the numerical model the mathematical problem needs to be properly posed. This requires the definition of the following:

• Model domain and boundary conditions
• Model layers and discretization
• Boundary conditions, internal sinks and sources
• Model parameterization
• Initial conditions and time stepping in transient simulations
• Model convergence

In developing (or reviewing) the numerical model setup, the modeller (or reviewer) should ask whether and how the numerical model boundaries and property distributions adhere to the conceptual model, and whether the model setup could potentially bias the model predictions in a way that is not intended.

This section provides an overview of the technical aspects of model setup (listed above) and provides guidance on their use in groundwater models for the natural resource industry. Where applicable, potential problems that may arise during model setup and their implications for model predictions are also discussed.

It should be recognized that the construction of a numerical model requires a solid understanding of the mathematical assumptions underlying the numerical modelling code and the various aspects of model setup. Although important, a discussion of the mathematical aspects of numerical modelling is beyond the scope of these guidelines. For more details on the numerical methods the reader is referred to standard groundwater modelling textbooks (e.g. Anderson & Woessner, 1992).

A comprehensive check list of specific information that should be considered for the numerical model setup is provided in Appendix D.

6.1 Model Domain and External Boundaries

The first step in model construction is the definition of a suitable model domain which should encompass the area for which groundwater flow (and potentially contaminant transport) is to be studied. The domain of the numerical model usually coincides with the domain of the conceptual model (but is sometimes smaller). The model domain encloses the geological volume of interest within which the mathematical simulation of groundwater flow and transport is computed.

The boundaries of the model domain are referred to as "external boundaries". The effect of these boundaries on heads and flows must then be conceptualized, and the best or most appropriate mathematical representation of this effect is selected for use in the model (also referred to as "boundary condition").

The delineation of the model domain is dependent on the selection of suitable external boundary conditions. It is generally preferable to use physical hydrological features which are known to control groundwater flow such as watershed divides, lakes, aquicludes etc. In some cases, the flow problem is simplified by defining a model domain which represents only a small portion of a larger, naturally bounded groundwater system (e.g. regional aquifer, basin, valley, island). When physical hydrologic features that can be used as boundary conditions are far from the area of interest, artificial boundaries are sometimes used. The external boundaries of the model domain must not be restrictive of the applied stress on the study site and must be compatible with modelling objectives (e.g. extent of dewatering or depressurization caused by mine pit excavation or pumping of wells). Therefore, whenever possible, the natural hydrogeological boundaries of the flow system and an appropriate model extent to capture the expected hydraulic stresses, should be used as the model domain boundaries of the numerical model (see section 6.3 for more details).

6.2 Model Layers & Discretization

A fundamental aspect of numerical models is the representation of the hydrogeological system with discrete elementary volumes. This discretization within the model domain is either represented by an orthogonal grid of model cells ("finite-difference" method) or a mesh of 3D elements ("finite-element" method) (Section 5.3). The spatial resolution or spatial accuracy of the model is limited by the size of the discrete volumes. In some cases the ability to represent a flow or transport process is also limited by cell or element size. Other considerations are the size of discrete features to be represented (e.g. fault zones, tunnels, dam faces, pumping wells, narrow geologic units, cutoff-walls, streams and ditches, etc.). A secondary consideration in cell or element spacing is the variability in aquifer properties - which is usually not known with greater accuracy than the typical model cell or element size - but becomes important in representing structural faults or boundary conditions in large regional models with large cell or element sizes. These guidelines provide examples of finite-difference models (with MODFLOW examples) and finite-element models (with FEFLOW examples).

The finite-difference grid has an orthogonal pattern and the model grid should be oriented to enclose the area of interest while minimizing the number of inactive (unused) model cells. If there is significant and uniform anisotropy in hydraulic conductivity across the site it may be advantageous to orient the model grid in this direction. However, most modelling codes allow the use of anisotropy at an angle different to the model grid. The finite difference grid can be refined by increasing the number or rows and columns abruptly or gradually (telescopic refinement) to obtain the desired refinement near features of interest (Figure 6-1).

The use of finite-element meshes provides more flexibility and can be more efficiently refined in many locations within the model domain (Figure 6-2). The finite-element mesh does not need to be oriented, and the mesh can be shaped to efficiently enclose a model domain. There are no inactive nodes in a finite-element mesh because the elements are fit exactly to the boundary. Anisotropy can be specified independently of a mesh shape.

Most models have a much larger horizontal extent than vertical extent, although in some situations the vertical extent may be as large as the horizontal extent (e.g. mine workings in very steep topography, models of vertical mine shaft dewatering). In most cases, the models are arranged such that the grid or finite-element mesh is used for model discretization in the horizontal or nearly horizontal dimension. However, for cross-sectional problems, the orientation of the grid or mesh can also be applied to a vertical plane to allow better vertical resolution. (e.g. idealized slope section of specified length represented in section-view by a finite-element mesh). Finite-difference grids result in rectangular volume elements, but the grid refinement can be specified only in one plane similar to finite-element mesh refinement in one plane.

Model layers are used for discretization in a direction perpendicular to the grid or mesh, usually in the vertical direction. In this sense, the finite-difference model, such as MODFLOW, is the same as the majority of commonly used finite-element model codes, such as FEFLOW. Both model types allow the use of layers that can be of uniform or variable thicknesses.

6.2.1 Cell or Element Size and Horizontal Discretization

Selecting the size of the cells or elements is an important step in numerical model design. The intended use of the model and the importance and size of the features being discretized affect both the evaluation of whether the model is discretized appropriately and whether important features are missing that would cause a systematic error or bias in the simulation results. Hydraulic properties and stresses are specified for each cell or element. Therefore, finer discretization of a model domain allows for finer definition of hydraulic properties and more accurate solution of stresses. In the past, the grid and mesh sizes were severely limited by the computational power of computers, but modern computers allow the modeller to use a very fine grid or mesh discretization that is often finer than required to solve a problem correctly. However, certain types of models, such as unsaturated flow or transport models, may require finer grid or mesh discretization.

The goal of discretization is to achieve an adequate model and prove it with sensitivity analysis, not to maximize model discretization to the maximum computer power available. Beyond certain discretization, the model solution does not change significantly, yet may be smoother in appearance or the water balance might be slightly better, but at a cost of longer computing times. At a minimum, the model cell or element size should be small enough to represent all the features of interest (e.g. mine tunnels, ditches, pumping wells, streams, shafts, small ponds, waste rock piles, etc.) as well as their properties and boundaries or sources and sinks. The model cell or element size should also be small enough to simulate applied stresses adequately. At the same time, the modeller should avoid introducing unnecessary cell discretization as this may unnecessarily complicate the modelling process. A larger, more complex model will take longer to run and may introduce more potential for error. In many instances, a simpler, smaller, faster model is preferred as it allows for more model runs during sensitivity and uncertainty analysis for a given modelling budget.

It is good modelling practice to start with a relatively coarse model grid and then subsequently refine the model discretization until the modelling results are stable and meet the required accuracy. The final model discretization should be discussed and justified in the modelling report. For projects, in which modelling results are sensitive to model discretization (e.g. flow and transport in heterogeneous media, such as bedrock with discrete high-permeability faults), sensitivity analyses on model discretization (e.g. doubling the grid or mesh refinement) should be conducted and discussed in the modelling report, with a focus on the results for the areas of increased grid refinement.

Local model refinement may also be required near important point sources or sinks such as pumping wells. The hydraulic head distribution is most complex near the well and the numerical solution is affected by the discretization of model cells. A finite-difference model does not simulate the hydraulic head gradient near a pumping well accurately, because the model extracts or injects water to the cell block volume which is usually much larger than the well diameter. The resulting head at the well is not a good approximation, but the heads away from the well are correct. If accuracy of the head near the well is not important to the problem, then the coarse grid is probably acceptable. But, if accuracy is needed near the well, then the finer grid would be necessary. Note that analytical equations are available to compute the drawdown in a pumping well simulated with a coarser grid. Finite-element models calculate the head in the pumping or injection well directly at the node, unless the well is placed between nodes. The hydraulic head is simulated more accurately at point sources, and sinks in finite-element codes.

Finite-difference models generally allow the widths of rows and columns to vary, which is called variable grid spacing. The use of variable grid spacing allows some flexibility to make cells smaller in some areas and coarser in other areas or to provide a smooth change of cell spacing). However, in some models there are advantages of using uniform, finely discretized grids, which provide the most accurate solution and simplify mapping of properties and boundaries. In finite-element models (e.g. FEFLOW), the mesh can be refined smoothly around points or lines or within specified areas, giving large flexibility in model design.

Figure 6-1: Examples of finite-difference grid discretization in MODFLOW Case Study 2.

Figure 6-2: Examples finite-element mesh discretization in FEFLOW.

6.2.2 Vertical Discretization

Vertical discretization may be required in a groundwater flow model to explicitly represent variation in aquifer properties with depth (e.g. the presence of a confining layer separating two aquifer units) or simply to provide better resolution of vertical gradients (even in a relatively uniform system). The degree of vertical discretization will depend on the modelling objectives. For example, a simple two or three layer model may be adequate for a regional flow model where horizontal groundwater flow dominates. In contrast, a detailed vertical discretization may be required to simulate inflow to an open pit in a bedrock setting with significant vertical heterogeneity.

Again, it is good practice to start with a relatively simple vertical discretization that covers only the essential features of the conceptual model (e.g. different hydrostratigraphic units). Sensitivity analyses can then be conducted to determine whether additional vertical discretization is required to adequately simulate the vertical aspect of groundwater flow. In any event, the model documentation should justify the vertical discretization used in the model.

Two approaches are commonly used to represent the discretization of a conceptual model in the vertical dimension: (i) uniform model layers (a rectilinear grid) and (ii) deformed model layers (see Anderson and Woessner, (1992) for more details). The deformed layers have layer boundaries usually following the surfaces of hydrostratigraphic units. Deformed model layers allow horizontal continuity to be maintained with fewer cells at the expense of introducing some error in the finite-difference and finite-element solution. The uniform layer approach has the advantages of simplicity and a stable solution, but may require more vertical discretization to adequately describe the vertical extent (thickness) of hydrostratigraphic units. It is also possible to create deformed layers for most units and to use regular (uniform) layers to subdivide selected hydrostratigraphic units achieving locally-fine vertical discretization (see Figure 6-3 from Case Study 2).

The model grid or mesh may also be used in profiling two-dimensional and three-dimensional models that have a rotated coordinate system (e.g. direction of gravity) to very accurately represent a portion of a site in a cross-section view. The limitations of such profile models are the inability to model radial flow to sources and sinks as well as the assumption of no flow boundaries as flow lines. However, such models are particularly useful in modelling a small section of approximately uniform slope or an elongated feature of interest.

Most three dimensional models, both finite-difference and finite-element use continuous layers in the vertical to model hydrostratigraphic units. Therefore, the problems with layer discretization are the same in these models (e.g. MODFLOW and FEFLOW programs). The following numerical issues should be considered when determining the vertical discretization in a layered numerical model:

• In deformed model layers only the dominant hydrostratigraphic units can be represented, and sub-units and heterogeneity is mapped by varying the unit properties locally (by zones). Where the layer surfaces are very irregular, the solution could be more difficult but the effect on numerical error is quite small (McDonald and Harbaugh, 1988). It is preferred to minimize layer irregularity and to keep layers horizontal or close to horizontal near planned mine excavations, where other inner boundaries will be added to the model (e.g. mine pit seepage face, mine tunnels, etc.).

Figure 6-3: Example of highly discretized finite-difference MODFLOW grid layers in vertical dimension with drain cells to simulate mine adit (Case Study 2).

• Intermediate layers should be used to separate larger layers representing hydrostratigraphic units of highly differing hydraulic conductivity in FEFLOW. These buffer layers increase model accuracy and have properties equal to one of the adjacent units.
• Pinching out of aquifers or confining beds in three dimensional models can be handled by changing the transmissivity of the layer or the hydraulic conductivity of the confining unit. In most cases, the hydrostratigraphic units are very simplified and modeled as either continuous layers or discrete units.
• In uniform and finely spaced model layers (Figure 6-4), the hydrostratigraphic units can be mapped with hydraulic conductivity distributions at any scale from regional to local. Small, nearly discrete features (e.g. mine shafts, adits, drains, and pit fill materials) are easily mapped onto a regular grid in three dimensions with the use of appropriate CAD or GIS software.

For density-dependent flow models, the grid size should be smaller than in flow-only models to achieve the optimum solution accuracy. In density-dependent models the horizontal size of each grid cells should be 0.38% of the model domain length, and layer thickness should be ideally 0.60% of model domain (Al-Kaktoumi et al, 2007). The solution accuracy may be acceptable at coarser grids and meshes, but it should be verified by the modeller.

Figure 6-4: Example of regularly spaced horizontal layers shown in vertical cross-sections through 3D model.

6.3 Boundary Conditions, Sinks and Sources

When evaluating the appropriateness of a groundwater flow model, the boundary conditions are critical, because they determine where the water enters and leaves the system. If the boundary conditions are inappropriate, the model will be a poor representation of the actual groundwater flow system. The modelling objectives and the magnitude of the stresses to be simulated also influence the selection of the appropriate boundary conditions. When groundwater systems are heavily stressed, the physical features that control the system can change in response to the stress. For example, a flow line boundary (described in Section 6.3.2.) determined from current groundwater levels may be a plausible boundary for simulating current groundwater flow but is often not appropriate for model predictions of future stresses such as pit dewatering. Any representation of these features must account for these potential changes, either by understanding the limitations of the simulation or by representing the physical feature as realistically as possible.

Boundary conditions are mathematical expressions of the state of the physical system that constrain the equations of the mathematical model (ASTM, D5609.1220124). In solving a groundwater flow problem, however, the boundary conditions are not simply mathematical constraints; they generally represent the sources and/or sinks of water within the system. Furthermore, their selection is critical to the development of an accurate model. Not only is the location of the boundaries important, but also their numerical or mathematical representation within the model, because many physical features, which are often hydrologic boundaries, can be mathematically represented in more than one way.

The determination of an appropriate mathematical representation of a boundary condition is dependent upon the objectives of the study (see Section 4). However, if the model is intended to forecast the response of the system to additional withdrawals that may affect the stage of the surface-water bodies, then a constant head is not appropriate and a more complex boundary is required. A model of a particular area developed for one study with a particular set of objectives may not necessarily be appropriate for another study in the same area with different objectives. All of these boundary condition parameters must be considered in evaluating the strengths and weaknesses of a groundwater flow model (USGS, 2004).

Water and contaminants can flow into and out of the numerical model through:

• External boundaries of the model domain.
• Internal sources and sinks (internal boundaries) within the model domain.

The flow system is usually a mix of specified head and specified flow boundaries. The use of flux boundaries only, without any specified head boundaries, should be avoided due to the non-uniqueness of the solution. Steady state problems require at least one specified head boundary to reference the hydraulic head solution.

6.3.1 Types of Flow Model Boundary Conditions

The boundaries can be generally grouped into two conceptual types: physical boundaries and hydraulic boundaries. Physical boundaries do not change in response to groundwater flow or stresses in the groundwater system, or at least the change is insignificantly small. Hydraulic boundaries depend on the groundwater flow system and may change as a response to groundwater flow.

• "Physical boundaries" are formed by the presence of a "large" body of surface water (see Section 6.3.3 for details), by the presence of a very low permeability hydrogeological unit (e.g. rock unit, fault filled with clay gouge, ice, or permafrost zone), by an artificial barrier in a small scale model, or by a natural drainage network (e.g. site located on top of a mountain surrounded by deep valleys).
• "Hydraulic boundaries" (artificial boundaries ) include groundwater divides and streamlines (flow lines) as well as distant head boundaries representing distant surface water bodies or arbitrary contours on a regional water table map. Hydraulic boundaries can be used to produce a steady-state flow field for calibration purposes, but these boundaries may not be acceptable for steady-state or transient predictive simulations of applied stresses. The use of a hydraulic boundary should be evaluated carefully to determine whether its use would cause unacceptable errors in the model.

All of the boundaries within the model domain, internal and external, can be represented with the same boundary conditions. There are many types of flow and transport model boundary conditions, where the hydraulic head, groundwater flux, or contaminant concentration is specified with some variations of constraints and interactions as well as variations in time and space.

6.3.2 No Flow Boundary (Streamline, Zero Flux)

A no-flow boundary does not allow any groundwater flow to occur across this boundary. In most numerical model programs, all the external surfaces and edges of the model domain are by default a no-flow boundary, if no other boundary type is specified. This means that the three-dimensional flow models do not "see" the ground surface and do not simulate any seepage of groundwater or runoff unless special boundaries are specified. The hydraulic head in any model node can rise or fall to any value, exceeding the vertical extent of the model if not constrained by boundaries. In most models, a no-flow boundary can also be specified anywhere in the model by using the specified flux boundary type (zero flux) or by using one of the hydraulic head boundaries with flux constraints specified as zero flux in or out of the boundary.

There are many common applications of no-flow boundaries:

• Groundwater divides are often chosen as no-flow boundaries, because they represent flow lines parallel to the boundary, so that there is no flow across the boundary. However, the locations of groundwater divides depend upon hydrological conditions and applied stresses, and can shift over time. Groundwater divides are dynamic boundaries and not physical boundaries of the flow system. Their representation as no-flow boundaries is justified only if the applied hydraulic stresses in the model domain have a negligible effect on the position of the boundary (e.g. the cone of depression of a water table or where depressurization does not reach the no-flow boundary during the applied stress such as mining excavations or pumping). See figure 6-5 for an example of where the stress reaches the boundary.
• Impermeable boundaries are natural geologic units or artificial materials which are regarded as effectively impermeable for modelling purposes, if the hydraulic conductivities of the adjacent materials differ by several orders of magnitude. In numerical models, a no-flow boundary can be represented with material properties, such as a very low hydraulic conductivity, without using hydraulic (head or flux) boundary conditions (head or flux). Mathematically, there is some flow in very low permeability materials and in models that do not allow zero permeability, but, effectively, such materials represent no-flow boundaries.
• In profile models of engineering structures, such as earth dams, the lateral boundaries are typically hydraulic streamlines, and the cross-section location in an elongated dam structure should be at least three times the depth of the structure. In other cases, the flow through the structure could be modeled in three dimensions.
• There are many complex engineering structures incorporating various cutoff walls, grout curtains, geomembranes and liners, frozen interiors, and other types. These structures could be considered as no-flow boundaries on a small scale, although most "impermeable" engineering materials are leaky to some extent and some may be very leaky. The choice of no-flow boundaries should be supported by engineering designs and field tests.
• The effect of density differences and transport of various fluids may in some cases affect hydraulic boundaries (e.g. a coastal aquifer saltwater interface or deep brine groundwater). For example, in a three dimensional flow model (freshwater flow only), the saltwater interface at the ocean shore is effectively a no-flow boundary with respect to water transfer because fresh and saline water does not mix easily in porous media. The use of a no-flow boundary assumes that there is no active salt water intrusion occurring. The sea water has tidal variation but if diurnal tides are not important for the project, the sea water body can be represented as a specified head of zero elevation along the shore line (see section 6.3.9. for details of specifying density-dependent boundaries).

Pre-mining End of excavation

Figure 6-5: Example of eternal hydraulic "no flow" boundaries affected by mine dewatering stress on hydraulic head maps from model simulation results.

6.3.3 Specified Head Boundary (Type 1 or Dirichlet condition)

At this boundary condition the hydraulic head is specified as a function of position and time. There are several variations of this boundary type:

• General specified-head boundary represents a hydraulic head which may change with time and in space. It is often used as a boundary where the specified hydraulic head is far away from the actual model cells or nodes where this condition exists. For example, these boundary conditions can be used to create a hydraulic gradient to a distant hydraulic head boundary (e.g. a distant aquifer boundary or distant surface water body) that is either constant or varies independently of the groundwater system (e.g. tidal fluctuations or river stages). This boundary allows use of a smaller model domain in a larger regional setting, while simulating hydraulic gradients appropriately between the local site being modeled and distant aquifer boundaries, which are located in or beyond the model domain. In most cases the specified head boundary is fixed in space at the boundary location, but it changes in time in transient models (e.g. mine pit excavation schedule, pit lake draining, or filling schedule).
• Constant-head boundary represents a hydraulic head which does not change with time and is not affected by the simulated groundwater system. It is often used to represent large surface water bodies which are not affected by stresses applied to the groundwater system in the model (e.g. large lakes or rivers, large streams in some cases, or the sea). Another use is to represent the observed hydraulic head (e.g. water table elevation) along a groundwater divide or even an arbitrary boundary sufficiently far away from the modeled site.

6.3.4 Specified Flux Boundary (Type 2 or Neumann conditions)

In this boundary condition the groundwater flux is specified across the boundary (node, line, surface) as a function of position and time. There are several variations of this boundary type (Figure 6-6):

• Constant flux boundary represents the simplest type of specified-flux boundary, where the flux across a given part of the boundary surface is uniform in space and constant with time. Practical applications in models vary, but it is commonly used to model known point or line source recharge (e.g. a leaking water management system component such as pipe or ditch) or discharge (e.g. the measured outflow from mine workings or boreholes).
• Specified flux boundary can be variable in time and is commonly used to represent a known time-variable point or line source recharges or discharges (e.g. a variable leakage from water management structures or injection/ or extraction rates from point sources, etc.). In 2D models, this boundary type is often used to specify the groundwater recharge flux into the top of the model (net precipitation). It is also possible to specify aerial recharge with a flux boundary, but this is an inefficient use of model boundaries (reduces options).
• General flux boundary represents a groundwater flux rate that is variable in space and time. This boundary condition is not commonly used or available.

6.3.5 Head Dependent Flux (Type 3, Cauchy or mixed boundary conditions)

A head dependent flux is a boundary type where a flux across this boundary surface adjusts in response to changes in the hydraulic head adjacent to the boundary. The flux is a specified function of the hydraulic head and varies during the problem solution just as the head varies. There are several types of hydraulic head-dependent flux boundaries:

• A very common use for this type of boundary is to simulate flow to and from rivers and drains (Figure 6-7). A separate river and drain boundary is available in some modelling programs (e.g. MODFLOW) to specifically represent groundwater and surface water interactions. The drain boundary behaves similarly to a seepage face boundary. In some modelling programs, this boundary is represented by a specified flux boundary with hydraulic head constraints (e.g. FEFLOW).
• This boundary may be used to represent flux through a leaky confining unit which is not explicitly represented in the model with material properties.
• Evapotranspiration is another example of a head dependent boundary. It is a flux from the water table which is often modeled as decreasing linearly with depth to water and becomes zero where the water table reaches some specified "extinction" depth. In most modelling programs, the evapotranspiration is implemented separately or as an option in the aerial flux recharge.

Figure 6-6a

Figure 6-6b

Figure 6-6: Examples of boundary conditions and internals sinks and sources involving specified flux: (a) recharge boundary zones in MODFLOW, (b) recharge flux and specified head boundaries in FEFLOW.

6.3.6 Seepage-face Boundary

A seepage face boundary is a mixed boundary condition (head-dependent flux) which is very often used in the resource industry's groundwater models. A surface of seepage is a boundary between the saturated flow field and the atmosphere, either by flux into the boundary from the model domain or by evaporation. A seepage face can only extract water from the model domain but cannot produce water into the model domain (one way flux). The seepage of groundwater occurs from the saturated portion of the model and the water table always intercepts the seepage face. The location of this type of boundary is specified, but the boundary may become active or inactive and change properties over time in some transient models (e.g. mine pit progressive excavation). The applications of this boundary type are for all slopes where groundwater discharge may occur (e.g. mine pits and tunnels, ditches, natural drain channels, dam slopes, waste rock pile slopes, ground slopes) - see Figure 6-7.

• In some modelling programs (e.g. FEFLOW), the seepage face boundary is a type of specified head-boundary with a constraint on boundary flux direction to the boundary (flow is only allowed toward the boundary). In other modelling programs, this boundary may be implemented as a separate boundary type. This type of boundary is used to simulate seepage from surfaces, linear features such as drain channels or tunnels, or points (e.g. drilled slope drains with screens or pumping wells with unspecified pumping rate).
• Other modelling programs may have a separate boundary type named Seepage Face Boundary , separate from all other boundary types (e.g. slope stability engineering models or some analytical 2D models) or the GUI may offer a separate boundary type that is internally solved using a combined head and constrained flux boundary.

6.3.7 Pumping Well (or Injection Well)

A pumping well is an internal source or sink and a mixed internal boundary type. Most three-dimensional modelling programs have this boundary type, and the use is exclusively for simulating pumping or injection wells and boreholes.

• This source/sink allows variable fluxes in time to be specified, often with hydraulic head constraints to turn off the flux if the head falls below a specified elevation. Well screens, well skin effects, and other options may be specified. Pumping wells use specified flow rates to extract (or inject) water, depending on the hydraulic head conditions. Pumping wells are used in three-dimensional or axisymmetric models to simulate radial flow toward the well.
• Two dimensional profile models cannot simulate radial flow and conventional pumping well source/sink "boundary types". Flux boundaries or drain nodes must be used, but the flow rates are not equivalent to real pumping rates calculated using three-dimensional or axisymmetric models.
• In some case it may be preferable to use a drain or a seepage face internal boundary at a point or line to simulate pumping or injection wells if the pumping rate is not known or is not important to the model solution. A seepage face or drain node allows the estimation of the maximum pumping rate possible at any point in time (e.g. for estimation of dewatering well requirements). This type of boundary allows for much faster transient model simulations because well drying and associated changes in time-stepping of the model are avoided.

Figure 6-7: Examples of mixed boundary type (specified head + flux constraints) to simulate an open mine pit seepage face and surface stream drains in a 3D transient flow model using FEFLOW.

6.3.8 Free-surface Boundary ("Water Table")

A free surface boundary is a moveable boundary where the hydraulic head is equal to the elevation of the boundary.

• The most common free-surface boundary is the water table. Some modelling programs allow the deformation of the model domain during model execution. Deformable free surfaces may increase model efficiency, but should be treated with caution in complex models to avoid unexpected effects and changes in model layers.
• Another example of a free surface boundary is the transition between freshwater and underlying seawater in a coastal aquifer. If diffusion and transient changes are not significant to the model objectives, the freshwater-saltwater transition zone can be treated as a sharp interface and a no-flow boundary of the fresh groundwater flow system.

6.3.9 Density-dependent Flow Boundaries

In numerical simulations where the water density varies in space and has a significant effect on groundwater flow, a density-dependent model must be used (e.g. FEFLOW, SEAWAT, SUTRA). The applications are in modelling of saltwater intrusion, waste injection into saline aquifers, heat storage in aquifers, brine disposal, geothermal flow, etc. (Anderson and Woessner, 1992). It is not common to see density-dependent flow models in resource extraction projects, but some coastal areas might experience salt water intrusion during excavation dewatering or pumping of aquifers, or there might be deep brines involved in the groundwater flow system at some locations (mostly in northern Canada and not likely in BC). The fresh water-saline water interface is used as an example in this section.

The geometry of the freshwater-saline interface must be estimated from observations, analytical solutions, or other numerical density-dependent flow model results. It is a type of a hydraulic boundary that may not affect short-duration transient models if the stresses do not reach the boundary (it then can be treated as no flow boundary). In long duration transient models, however, the fresh water-saline water boundary may also change with time. In some cases, the model can simulate the formation or dissipation of a freshwater "lens" or zone above the saline water. Saline water may also be discharged into fresh water and then sink down due to density differences. The modelling objectives determine how these boundaries are specified.

The hydraulic head boundaries must be expressed as equivalent freshwater head values. The head nodes along the sea shore, for example, are assigned zero head values (or some exact mean sea level or transient sea level schedule). In three dimensions, or in two dimensional profile models, the equivalent freshwater hydraulic head varies with depth below sea level (or some reference elevation) as a product of density ratio between fresh and saline water and depth. This means that at appropriate model layers or slices or nodes, the specified head boundary values vary with depth (see example of an axisymmetric model with this boundary type in Figure 6-8).

Note that, as with transport modelling (see Section 9), there are additional parameters required to simulate density-dependent flow: representative mass for each water type (represented with chloride concentration or other), diffusivity and dispersivity, both of which are generally low. The reference mass is assigned for the concentration of a reference parameter (e.g. chloride) in fresh water originating from freshwater recharge flux source. The recharge flux boundary requires a constant "freshwater" concentration boundary to be specified as source of freshwater.

Corresponding to the specified head boundary of the fresh water-saline water interface is a constant concentration boundary with mass flux constraints to allow any types of water to exit the model (seep to shores) in the top model layer or slice, but allows the saline water flux to enter in other model layers or slices (saline water intrusion). For more details on density-dependent flow modelling the reader is referred to Anderson and Woessner, 1992; USGS, 2002, WASY, 2007).

Figure 6-8: Examples of boundary conditions in axisymmetric model (2D cross-section) in density-dependent flow model.

6.3.10 Uncertainties of Boundary Conditions

The choice and setting of boundary conditions is a critical step in model design. Errors in boundary conditions may cause serious errors in model predictions. The uncertainty associated with boundary conditions can be described in four broad terms:

• Conceptual model representation of boundaries, sources, and sinks
• Implementation of boundaries in a numerical model
• Numerical solution (e.g. a water table)
• Stress dependency of boundaries (see next section)

Conceptual model representation

The physical or hydraulic boundaries must be represented appropriately and agree with the conceptual model in order to achieve good model results (see Section 4). The boundaries must be explained by the real hydrogeological system and supported by some observations (e.g. water levels) or by reasonable inferences based on physical process understanding (e.g. catchment drainage network and ground topography vs. groundwater divides). The known properties of geologic materials or artificial materials (e.g. waste rock dumps) should be reviewed in the conceptual model part of the modelling process, and the boundary conditions must take the whole hydrogeological system into consideration, including surface water and any planned future changes to the surface water drainage pattern.

Boundary condition implementation

Once the appropriate types of model boundaries, sources, and sinks are selected, these boundaries must be appropriately implemented or coded in the modelling model to achieve the intended behaviour as specified in the conceptual model. Common errors include:

• Data entry errors (e.g. errors in head or flux values or errors in the direction of flux)
• Wrong assumptions are used that are not supported by the conceptual model (e.g. specified head node or cell is an unlimited source of water)
• Oversimplification of boundary geometries relative to model objectives
• Too excessive complexity in boundaries that is not supported by observations (e.g. complex distributions of recharge fluxes to achieve better calibration of the heads)

The modelling program must correctly solve the flow and transport simulation and must properly use the boundary conditions. Numerical instabilities and lack of convergence may still occur in some cases in localized areas of the model domain, even after model convergence of the specified criteria for head and flux. In rare cases, the modelling files or the GUI may be corrupted and unintended behaviour of boundary conditions observed (or missed). Modern modelling GUIs present the hydrogeologist with many options and settings, which should be understood before using, and/or investigated numerically in a sensitivity analysis. The numerical behaviour of boundary conditions must be evaluated in terms of reasonable and expected behaviour at all steps in the numerical simulation process by the modeller.

Groundwater recharge is a very uncertain flux boundary. Typically, there is little data for total recharge at a site, much less for spatial distribution of recharge. However, recharge rates and zonation is usually defined through model calibration to observed heads, given some known aquifer properties. The uncertainty of aquifer properties will directly affect the uncertainty of the "calibrated" recharge flux, and the model solution is usually non-unique. For example, in a steady-state model, any combination of recharge and hydraulic conductivity can result in an identical head solution if the ratio of recharge to hydraulic conductivity (R/K) is the same. The difference between different combinations will be observed only in the model water balance.

Numerical solution of a water table boundary

The water table is a special boundary which is almost always solved numerically and not pre-specified. The models may use unsaturated/saturated conditions and represent the water table boundary as the surface of the zero pressure head. However, unsaturated parameters are difficult to estimate and the numerical solution has many complexities and is numerically demanding. Many three-dimensional models (e.g. MODFLOW) use Dupuit assumptions to approximate flow in the top (unconfined) model layer.

6.3.11 Stress Dependent Problems

The model boundary conditions must be examined if the model objectives include modelling of applied stresses to the groundwater system. The model representation of a system boundary may be a function of the nature and magnitude of stress applied to the system during model simulation. If the boundary conditions are stress dependent, the model cannot be considered a general, all-purpose tool for investigating any stress on the system, because it will give valid results only when the stresses do not impact the boundary. The study of a new stress on the same model may require the reformulation of the representation of model boundaries and sensitivity tests on the model boundary representation.

Stress-dependency is a primary concern wherever the model boundaries differ from the natural (physical) system boundaries. For example, the modeller may extend the model boundaries beyond the groundwater divide (i.e. a natural boundary) to simulate drawdown due to long-term pumping. Stress dependency of external boundaries of a model domain is usually evaluated with a larger model (larger domain extent) before a detailed local model is constructed. Each internal boundary can also be checked through a sensitivity analysis (adding, removing, or changing of boundary conditions and observing effects on the model simulation from the applied stresses).

Much of the sensitivity analysis to boundary effects in the model is done by the modeller during model development and calibration, and most of the problems are corrected without documentation in the modelling report. In cases where boundary interactions are observed in the predictive results, or if the model results are unexpected, a full sensitivity analysis should be presented on the boundary effects.

Boundary conditions strongly influence modelling results and the modeller (and reviewer) should critically evaluate the potential for errors in defining boundary conditions. The following examples demonstrate common errors:

No-flow boundaries affect stress calculation:

• No flow boundaries along the model's edge may also cause mounding of recharge to levels above the ground surface in unconfined aquifer layers of the saturated flow model (see Figure 6-9a).
• No-flow boundaries of the model domain affect the model's response to applied stresses at the site of interest. A hydraulic no-flow boundary is usually positioned far enough from the area of interest that whatever hydraulic stresses are simulated in the area of interest will not reach the distant hydrological boundary. A significant change in the hydraulic head near the no-flow boundary suggests that the stress has reached the boundary and the solution is beginning to be affected by the boundary presence. Hydraulic head contours of a cone of depression caused by dewatering will be "bent" by the no-flow boundary on the head contour maps. An example of a drawdown caused by dewatering near an excavation and the effect of no-flow boundaries is shown in Figure 6-9b). In most models, the no flow boundaries are external boundaries only, therefore, the effect of no-flow boundaries, if present, will depend on the aquifer properties, distance to the area of interest where a stress is applied, magnitude and duration of applied stress, and other factors.

Specified head boundaries produce too much water:

• A more serious problem may arise from inappropriate use of specified head boundaries because the model may gain an unintended large flux of water from unconstrained specified head boundaries. A single inappropriately placed and unconstrained head node in the model can unintentionally (and erroneously) produce large quantities of water and completely change the water balance of the model solution. In this situation, the specified head boundary will artificially constrain the flow system by providing or accepting fluxes as a result of the nearby applied stresses.
• A classic example would be a small stream or pond which produces a large quantity of water into the model, much larger than is possible from its small observed stream flow or pond water balance. It is difficult to detect these errors during initial model review, but the water balance of the flow system should be examined for groundwater sources and sinks and any unusually large fluxes (see Figure 6-9c).

Specified head boundaries show unrealistic behavior during some stresses:

• A specified head boundary may function appropriately in natural conditions, but have unrealistic behavior under stress conditions. For example, a small to medium-sized stream, may function as a specified head boundary if the stress does not induce flow to or from the stream of sufficient magnitude to significantly affect the stream stage. If, however, the stress is so large as to cause a part of the stream to dry up, then the stream can no longer be treated as a specified head boundary. The stream may need to be modeled as a flux-dependent head boundary.

Specified flux boundary lacks head-dependency constraints:

• Specified flux boundaries without hydraulic head dependency assume the flux to or from the model is independent of the stress in the model. For example, applied recharge flux will not normally stop the flux when the water table rises to the ground surface leading to artificial mounding of the head above the ground surface (i.e. flooding) of the model. This behavior may not be noticed if there are many nearby drains which accept the recharged water, but the water balance of the model will be affected and may be in error. In particular situations, a head-dependent flux may be more locally appropriate in sensitive model areas.
• Spatially distributed recharge fluxes should be justified based on observations, soil properties, surface water hydrology, and climate. Highly variable recharge fluxes over the site of interest, specified to help model calibration, should be questioned. Any model can be calibrated to the head conditions using variable recharge distribution, but such recharge variability must be physically justified.

In summary, the model boundary conditions must be examined if the model objectives include modelling of applied stresses to the groundwater system, especially when hydraulic or artificial boundaries are used. Great care must be taken in implementation of boundary conditions in numerical models and the model results, including water balance, must be examined.

Figure 6-9: Examples of potential problems with specified head boundaries in a homogeneous aquifer.

6.4 Model Parameterization

6.4.1 Concepts

An important aspect of representing the hydrogeological conceptual model is the choice of the hydraulic properties assigned to the model domain. In the saturated flow model, the properties include hydraulic conductivity and or transmissivity, unconfined and confined aquifer storage properties (specific yield and specific storage), and aquifer compressibility. In the unsaturated flow models, there are more properties specifying the function of hydraulic conductivity with saturation (i.e., soil water characteristic curves). Transport models will require additional properties to describe the interaction of contaminants with porous media (see Section 9).

In every numerical flow model, the properties are assigned to model cells or elements and are usually assigned to layers of cells or elements. Each model cell or element can have different properties. Given the uncertainty of knowledge of the distribution of hydraulic properties from available data, the modeller must choose the most appropriate representation of these properties. The following types of spatial distributions of properties are commonly used in flow models:

• Uniform distribution is based on a single representative value of a parameter in a large unit or in many discrete sub-units each with uniform distribution, (i.e., one value).
• Interpolated (spatially variable) distribution is based on an assumed or observed variation of the parameter in space.
• Anisotropy , specified for hydraulic conductivity, is based on geological information and sometimes on hydraulic testing results. It is often used to deform the flow field to help with model calibration.

6.4.2 Uniform Distributions

The mean value of hydraulic testing results is usually used to represent one homogeneous hydrostratigraphic unit or each sub-unit within a larger unit (Figure 6-10). The method of averaging hydraulic conductivities is not a simple choice, and it depends on the type of statistical distribution of data, flow direction (e.g. vertical across many hydrogeological units or horizontal within one hydrogeological unit), and the nature of heterogeneity of porous media (or equivalent porous media for fractured bedrock). There are three commonly used averaging methods to produce a mean value, which increases in value in this sequence: arithmetic mean, geometric mean, and harmonic mean. In most models, the hydraulic conductivity is adjusted during model calibration, but the field test data usually provides reasonable bounds for this parameter.

Most hydrogeological studies and modelling results support the use of the geometric mean of K, for bedrock aquifers assuming a good site-wide connectivity of fractures, although some hydrogeologists prefer the harmonic mean for fractured rocks in a poorly connected fracture network. The geometric mean is used where the hydraulic conductivity distribution is log-normal (or normal), and the distribution is spatially random (i.e. hydraulic tests are independent). The harmonic mean is used to average the hydraulic conductivity values for flow that is perpendicular to different hydrostratigraphic units (e.g. flow across a low permeability fault or vertical recharge across a confining unit). Overall, the choice of the representative value or a mean value of model parameters has a large influence on model results and predictions and should be reviewed carefully.

Figure 6-10a

Figure 6-10b

Figure 6-10c

Figure 6-10: Example of vertical discretization of hydraulic conductivity of hydrostratigraphic units: (a) MODFLOW model, (b) FEFLOW models of mine site with deformed, and (c) uniform layers.

6.4.3 Interpolation of Spatially Varying Properties

Interpolation of heterogeneous property distributions should only be done if there is a sufficient (large) data set and strong geological data supporting such spatial distributions, and if the interpolated property will not be adjusted during model calibration. For example, a tailings pile could be tested on a regular grid in many test holes and the resulting hydraulic conductivity interpolated. Generally, the accepted practice is to use the simplest interpolation method that is consistent with the known data. However, it is possible to construct different modelling scenarios with different parameter distributions resulting from various interpolation methods fitted to the data. The spatial variation can be calculated using various interpolation methods (e.g. distance weighed, linear, polynomial, or radial basis functions) or using more rigorous geo-statistical methods (e.g. kriging). This is rarely done except for research purposes. The simplest interpolation methods (e.g. linear trend over the study site) are preferred.

An interpolated hydraulic conductivity distribution is very impractical to calibrate, because the whole distribution must be changed in each calibration step. Therefore, a uniform distribution of K values within each unit or sub-unit is preferred for calibration of K values.

6.4.4 Anisotropy

Anisotropy is almost always specified over the whole hydrostratigraphic unit or locally within a specified model volume. There is almost never enough data to smoothly interpolate a spatial distribution of anisotropy. The orientation of the principal axes of anisotropy has a large effect on flow simulations, especially in steep terrain. It is important to recognize that there are different principal axis configurations for anisotropy:

• Vertical and horizontal in Cartesian coordinates - This anisotropy is aligned with geographic directions and the direction of gravity (Kx, Ky, Kz). It could represent horizontal bedding or fracturing of geological units that have been eroded or intruded, but where the anisotropy within the unit remains relative to horizontal or vertical axes.
• Parallel to slope of model layers - This is useful for deformed sedimentary geologic units or various types of sheet fracturing due to unloading of stress that occurs parallel to existing slopes and irregularly shaped units.
• Anisotropy specified in three dimensions using angles for each axis - This anisotropy direction is independent of Cartesian coordinates or model layer orientations and may be required in bedrock settings with significant tectonic movement.

It is common in models of alluvial aquifers and other sedimentary rock aquifers to specify vertical anisotropy of some assumed (or typical published) value. There should be a good discussion of the choice of anisotropy ratios and sensitivity analysis of this input. In steep topography where large gradients exist, anisotropy perpendicular to slope allows groundwater to flow easier in a downslope direction, while limiting water recharge and discharge perpendicular to the slope. The three-dimensional geometry of the model will determine how much of an effect anisotropy has and the sensitivity analysis will demonstrate whether the model results depend on, or are strongly influenced, by anisotropy assumptions.

6.4.5 Aquifer Heterogeneity

An important consideration in model parameterization is the representation of the natural variability of the aquifer system, i.e. the degree of heterogeneity and anisotropy of hydraulic properties (see discussion in Section 4). The degree of heterogeneity that should be incorporated in the model will depend on the complexity of the site and the modelling objectives but also on the availability of field data. At study sites with very sparse data, the modeller has no choice but to use homogeneous and isotropic assumptions. At sites where heterogeneity and/or anisotropy are well documented and important for the modelling objectives, this it should be explicitly included in the numerical model.

In most groundwater models for natural resource projects, aquifer heterogeneity is only considered at the domain scale. In other words, a uniform distribution of K values is assumed for each hydrostratigraphic units (or sub-units) (e.g. using an average K value). Smaller-scale heterogeneity, i.e. variations in K at the scale of an individual hydrostratigraphic unit (or subunit) can be modeled using geostatistical methods (see above), but requires detailed field information and is usually not considered for the large-scale models used in natural resource projects.

6.4.6 Aquifer Storage Properties for Transient Simulations

In transient simulations of groundwater flow, the storage properties of hydrostratigraphic units are used to determine the volume of released or stored water at each simulation time step. Unconfined storage is related to the release of water as the water table lowers (dewatering of the aquifer material); thus, it occurs only along the top boundary of the saturated flow system. Confined storage is related to the release of water as the head drops, because of the expansion of the water itself as the pressure changes. It is also related to changes in the solid skeleton of the aquifer (no dewatering occurs). In simulating changes in storage for transient systems, it is important that the unconfined storage occurs only at the top boundary (or top active layer), even if the water-table aquifer is divided into many layers. Most three-dimensional models (MODFLOW and FEFLOW) automatically control which storage coefficient is used based on the layer geometries and heads, thus ensuring the proper coefficient is used. Some model programs may require the modeller to specify the coefficient for each cell.

Storage properties are based on large scale pumping test results with at least one observation well. Not all natural resource sites have pumping tests completed, and in those cases, the storage parameters are taken from published properties for similar geologic formations or at nearby sites. These parameters can be adjusted through transient calibration, but should be included in sensitivity analyses.

6.4.7 Evaluation of Specified Properties

The selection of adequate parameter values and spatial discretization of hydraulic parameters for a groundwater flow simulation requires experience in hydrogeology and numerical modelling. The criteria include: the continuity of deposits, the type of depositional environment, the type of fractured rock, and the type of test data available, among others. As always, the objectives of the study also determine which features must be represented in the model and the level of detail required to adequately represent their effect on the flow system.

6.5 Transient Simulations

Transient models simulate time-dependent problems such as the impact of stresses over time. A transient simulation must begin with representative initial conditions and ends at a specified elapsed model time, which depends on the modelling objectives and time scale of the problem. Time is divided into time steps, and hydraulic head is computed at the end of each time step. Transient models are inherently more complicated than steady state problems, take longer computation time and create much larger data storage demands.

The special demands of transient simulations are as follows:

• The storage properties of the hydrogeological units must be specified.
• Initial conditions of hydraulic head distribution must be defined and must represent the site conditions.
• Boundary conditions (external and internal) must be adjusted to transient conditions (e.g. seasonal recharge flux rate, seasonal river stage variation, sea tides, changes in depth of mine pit or tunnel excavation, changes in mine site surficial materials, diversions of natural streams).
• Boundary conditions must be checked for interactions with hydraulic stresses which propagate over time especially at later times of the simulation.

6.5.1 Initial Conditions in Transient Simulations

Accurate definition of initial conditions for transient groundwater models is an essential part of conceptualizing and modelling transient groundwater flow (ASTM, D5610.1220124). Initial conditions represent the heads at the beginning of a transient simulation, and these conditions serve as a "boundary condition" in time for the transient groundwater model solution. Initial conditions for a flow system are usually represented by the hydraulic head distribution throughout the model domain (e.g. pre-mining conditions, average aquifer conditions before pumping). Initial fluxes may be also specified (e.g. pre-mining recharge rates). As the model simulation progresses, the new calculated changes in the hydraulic head through time will be relative to these initial heads.

Initial conditions are used only in transient simulations and are different from starting heads (or the initial guess) in steady-state solutions. In steady-state solutions, the starting heads affect the efficiency of the numerical solution, but the final solution is not dependent upon them; i.e. the same result should be observed independent of starting heads. In transient solutions, however, the initial conditions are the heads from which the model calculates changes in the system due to the stresses applied. Thus, the response of the system is directly related to the initial conditions used in the simulation.

There are two commonly used methods for defining initial conditions for transient simulations of hydraulic stresses: defining or simulating the steady state conditions or simulating long-term transient conditions prior to the modeled stresses in the predictive model.

Defining steady-state initial conditions is the preferred method if the hydrogeological system is in an approximately steady state prior to expected stresses. Any natural ground conditions without recent large disturbances to surface water drainage or extreme climatic change are typically assumed to be in steady state condition on an inter-annual time scale. Alternatively, a steady-state can be demonstrated with measurements. The numerical model can be solved in steady-state mode, and the steady-state solution imported into the transient model as part of the initial conditions. A steady-state model should be calibrated to the steady-state conditions.

Defining transient initial conditions is required if the hydrogeological system is not in steady-state equilibrium. This condition could occur (i) if the simulation is intended to simulate seasonal or other cyclic conditions where the system is never at a steady state (ii) in instances where there is a period of unknown stress that cannot be reproduced accurately, or (iii) when it is not feasible to simulate the entire period of record from a time of steady state. Under these conditions, it is important the initial conditions used do not bias the results for the period of interest.

The use of model-generated head values ensures the initial conditions and the model setup are consistent. In contrast, the use of field observations (say groundwater level time trends) for initial conditions can be problematic, in particular if the assumed initial head distribution will be inconsistent with the model setup (model setup precedes model calibration so there might be unanticipated adjustments in hydraulic heads once the model is run in transient mode).

In some cases, where the transient initial conditions are not representative, the model solution may be affected and the modeled response to stress may include a numerical model solution adjustment from initial time (at initial given conditions) to the model's boundary conditions - and not only to the applied stresses. The following situations are typically susceptible to this type of error:

• Slow and small expected response to applied stresses (e.g. aquifer is of low permeability, or stresses are small in magnitude or both).
• Highly time-variable hydrogeological system (e.g. highly-seasonal recharge and/or high variations in hydraulic heads).
• Changes in the hydrogeological system which are not accounted for in the initial conditions of the actual site, which may have occurred since the steady-state data was collected (e.g. river diversion or control, strong climatic change trends, or changes in groundwater use near the study site).

In other cases, the model is not significantly affected by initial conditions due to the high permeability of the aquifer and rapid equilibration to any stresses, whether natural or caused by natural resource extraction activities.

The determination of the importance and duration of erroneous effects or imperfect initial conditions can be accomplished by testing the effect of different initial conditions on the model under study. This test is accomplished by simulating the same system with the stresses and different initial conditions. When the simulations for all the different initial conditions produce the same result, then one can assume the influence of the inaccurate initial conditions is negligible for all the following time periods. There are published time constant formulas for estimating the aquifer response to transient stresses (Domenico and Schwartz, 1998), but it is recommended to model this response numerically for natural transient conditions of the site and to present graphically the results in support for the choice of initial conditions.

6.5.2 Initial Conditions in Transient Density-dependent Flow Simulations

There are two types of initial conditions: concentrations (densities) and equivalent freshwater heads. At each model node, the initial water density must be specified (sometimes specified as reference chloride concentration or another parameter concentration). At a maximum reference concentration, the density of water is assigned equal to that of saline water (as specified by the density ratio). At minimum reference concentration (background concentration), the water density is that of fresh water. Intermediate values are calculated numerically between these bounds.

6.5.3 Time Stepping

The intended use of the model is an important factor in choosing the length of stress periods and time steps. Many time steps are required to simulate a complex distribution of the head over time. Also, the size of the time steps has an impact on the accuracy of a model, particularly at the early stages of each applied stress. This is similar to the need for many cells to represent the spatial distribution of hydraulic head. It is important to incorporate enough time steps to allow the temporal complexity of the head distribution to be simulated (USGS, 2004).

Modelling programs such as MODFLOW and FEFLOW use stress periods to represent stress data as a function of time. A stress period represents a group of one or more time steps in which stress input data (e.g. pumping rates) are constant. A new stress period must start whenever it becomes necessary to change stress input data. If stress periods are too long, important dynamics of the stresses may be left out or poorly represented. For example, the Well Package of MODFLOW allows pumping rates for wells to change every stress period where the pumping within a stress period is constant (stepped function). In FEFLOW, the pumping rate can vary linearly or non-linearly over time, but the model solution is much slower due to small time steps used to vary the applied pumping stress incrementally.

There is a balance between time stepping accuracy and modelling time (and budget). For example, open pit excavation could be done progressively with very fine increments and appropriate geometry of pit depth and shape to achieve a precise solution. However, the simulation time would be very long and there is uncertainty about the actual excavation schedule which normally varies by many months, and for operational reasons, so does the pit geometry. The models are usually done by pit "stages". Each stage is introduced at some discrete time during model simulation and the transient model shows how the groundwater system responds to each stage of excavation. A stepped approach of large applied stress is a simplified and reasonable representation of that type of boundary condition.

There is a need for performing a sensitivity analysis of the transient simulation results to the time stepping scheme of the model and the time-discretization of the applied stresses. Most modellers make several trial model runs with different time steps to check for effects on the final solution, but a formal sensitivity analysis may not be done. If a model is used to analyze the average response of a system over many years, then the applied stress might be represented as yearly averages using yearly stress periods. If a monthly or seasonal analysis is required, then the stress periods should also be small enough (say monthly or even biweekly) to represent that temporal variation in the recharge, the other boundary conditions, and the applied stresses.

There are published criteria for choosing the maximum time-step duration, which depends on the square of the model or element cell area and on "constants" such as the aquifer storativity and transmissivity. Most modelling programs have either automatic time stepping control, based on such criteria, or allow user-specified time stepping progression (1.2 to 1.5 is considered appropriate). It is the responsibility of the modeller to check the sensitivity of the modelling results to the time-stepping settings.

6.6 Model Convergence

Groundwater numerical models calculate the solution of large sets of simultaneous algebraic equations using sophisticated matrix solution techniques ("solvers "). Most of the solution techniques are iterative, where the solution is obtained through successive approximation, until a "good" solution has been obtained (the model converged on a solution). The criterion used in most iterative solution techniques is called the "head change criterion". When the maximum absolute value of head change from all nodes during an iteration is less than or equal to the selected head change criterion, then iteration stops.

6.6.1 Accuracy of the Matrix Solution

When evaluating groundwater flow model results, the reviewer depends on model documentation and the reported head change criterion at model solution convergence. One means of evaluating whether an appropriate head change criterion was used is to examine the global mass balance for the model. If the error in the mass balance (e.g. total inflow minus total outflow divided by one half the sum of the inflow and outflow) over the entire model domain is small, usually less than 0.5 percent, then the head change criterion is assumed to have been sufficient. If the error in the mass balance calculations is significant, then the model solution was not good enough.

Even if the head change criterion is met and the global mass balance error is small, the model solution may not be appropriate for the system under investigation. Two potential reasons are that some models can either be mathematically non-unique or very nonlinear. The mathematically non-unique problem is usually a poorly posed problem, where a model has only specified-flow boundary conditions and no other boundary condition that specifies a head or datum (such as, constant head, river stage, general head boundary, etc.). In this type of problem, there is a family of solutions all with the same gradients, but different absolute heads. The matrix solution technique may not converge or it may converge to one of the infinite number of possible solutions.

The accuracy of the matrix solution is usually not an issue with groundwater models that meet the head change criterion and have small mass balance errors. It is important when using nonlinear models. In nonlinear problems, the solution affects the coefficients of the matrix being solved; thus, the solution affects the problem being solved. As a result, the manner in which the iterative solution technique approaches a solution can affect the final solution. The rate of convergence and the method of making cells inactive is usually considered by the modeller, but it is not easy to determine if the solution is correct.

6.6.2 Methods Used to Improve Solutions

The matrix solution can be improved by lowering the head change criterion, adjusting iteration parameters (if the solution techniques use iteration parameters), using different starting heads for steady-state simulations, or using a different solution technique. Some models do not converge smoothly, and modellers use non-standard methods to obtain a model solution. As long as the non-standard method does not violate any important hydrological processes, they are usually transparent to the final solution and are appropriate. However, these non-standard techniques should be evaluated to determine whether they cause potential errors to be introduced into the model solution.

Examples of non-standard modelling methods to improve a solution include:

• Saving of intermediate solutions that have not yet converged and changing matrix solution parameters when restarting the model.
• Making a nonlinear water-table simulation linear by fixing the saturated thickness of the model to avoid cell re-wetting and drying problems as well as head oscillations.
• Obtaining a steady-state solution by using storage to slow convergence and dampen the approach to the solution through simulating a long transient time period.

Review Questions

1. The most important reasons for model mesh or grid discretization are:
   1. To allow appropriate representation of the most important features in a groundwater system and distributions of hydraulic properties.
   2. To provide the most detailed model output which can be achieved given available computing power.
   3. To decrease the model run time.
   4. To provide fine model grid or mesh density to guarantee precise model results.
   5. A and D.
   6. All of the above.
2. Definition of boundary conditions is one of the most important modelling steps because boundaries determine where water enters and leaves the system and thus, proper selection is critical to the development of an accurate model.
   1. True.
   2. False.
3. Which of the following boundary types can best represent a large surface water body?
   1. Specified flux boundary.
   2. Specified head boundary.
   3. Head dependent flux boundary.
   4. No flow boundary.
   5. None of the above.
   6. All of the above.
4. The main differences between steady state and transient simulations are:
   1. In transient simulations, the storage properties can be used to determine the volume of released or stored water at each simulation time step.
   2. Steady-state models do not use storage properties.
   3. Transient models are more complicated because there is one more dimension to the model (time).
   4. Accurate definition of initial conditions is an essential part of conceptualizing and modelling transient groundwater flow.
   5. All of the above.
5. Which of the following does NOT affect the numerical accuracy of a converged solution in numerical model (assuming a good conceptual model and calibration)?
   1. Water mass balance error.
   2. Precise presentation of model results in documentation.
   3. Time stepping in transient model.
   4. Model domain discretization.
   5. Head convergence criteria.
   6. Poor value of a calibration average error measure such as nRMSE.
   7. B and F.
   8. All of the above.

---

Proceed to Section 7: Model Calibration & Verification