Mathematical Model Selection

5.1 Overview

Once a conceptual model has been developed (Section 4) the groundwater modeller has to select an appropriate mathematical model to represent the conceptual model. This process encompasses (i) the selection of analytical versus numerical models; (ii) the selection of spatial and temporal dimensionality for numerical models; and (iii) selection of a computer code or analytical solution approach.

The selection of a mathematical model involves consideration of the following (see Figure 5-1):

• Modelling objectives
• Conceptual model
• Modelling complexity
• Modelling tools available
• Modelling constraints.

The most important factors to consider in mathematical model selection are modelling objectives and the conceptual model. The modelling objectives and site conceptual model are important because they determine what processes should be simulated and what level of accuracy is required. The required accuracy determines the model complexity which in turn influences the selection of the mathematical model, and the extent of data collection used to drive the model (Figure 5-1). For example, an analytical solution may be adequate for order-of-magnitude estimate of steady-state inflow to an open pit. In contrast, a three-dimensional numerical model may be required if the modelling objective calls for an estimate of pit inflow to an accuracy of +/- 25%.

The conceptual model provides a framework for selecting the appropriate level of modelling complexity and ultimately mathematical model selection by defining the relevant processes that must be included (Figure 5-1).

Note that the process of mathematical model selection is closely related to the development of a conceptual model (Section 4) but usually requires a greater level simplification due to mathematical constraints. These guidelines therefore distinguish mathematical model selection from the more qualitative model conceptualization.

This section introduces the key definitions and concepts and provides guidance on appropriate model selection for natural resource projects. This section includes a discussion of the types of models commonly used, their capabilities and limitations, and their applicability to various physical conditions and applications.

Figure 5-1: Factors influencing selection of a mathematical model.

5.2 Scope of Mathematical Model

The first step in model selection is the definition of an appropriate scope of modelling. This scoping step will drive the decisions regarding analytical versus numerical modelling, dimensionality of the model, and selection of an appropriate code or solution to address the modelling objectives. The following questions have to be evaluated during model scoping:

• Should the groundwater flow problem be simulated in three dimensions or can it be simplified to a 2-dimensional problem, an axi-symmetric (radial) problem or even a one-dimensional problem? Can analytical approaches be used (i.e.1D or uniform radial flow)? Is vertical flow important to model?
• Are conditions over time important to model? Should the groundwater flow problem be simulated using a transient model or can it be approximated by a steady-state model? If so, what steady-state conditions should be modeled (e.g. high flow, low flow or average flow, or a representative range)?
• If the setting is bedrock, would an equivalent porous medium approach be appropriate and can this approach be justified, based on site characteristics? If not, would other approaches like discrete fracture network or dual porosity model be required and what data from the site is needed to support these models?
• Can flow in the unsaturated zone (or in perched zones) be ignored or is a variably saturated flow model required?
• Should contaminant transport be simulated using a contaminant transport model or is the use of particle tracking adequate?
• Should geochemical aspects of contaminant transport be simulated or is the use of a conservative solute transport model adequate?

The above questions have a strong influence on the degree of complexity of the mathematical model and often the type of model (i.e. code) required for the modelling study.

The answer to these questions depends on the modelling objectives and the required model accuracy. The selected scope of modelling depends on the conceptual model and the available field data on which the conceptual model is based.

The following sections provide guidance on the appropriate scope of groundwater modelling to assess impacts on natural resource projects.

5.2.1 Model Dimension

Spatially, groundwater flow systems are three-dimensional; that is, they can have a significant flow component along each of the three orthogonal axes, x , y , and z .

Using the principle of parsimony or simplicity (see Section 4) the modeller should seek to simplify the flow problem as much as possible. This includes an assessment of whether the flow problem can be reduced to a 2D or even 1D problem (if flow component is not significant in 1 or 2 of the other dimensions). The following sections discuss the criteria that should be considered in determining the appropriate dimension(s) to be used in the groundwater flow (and transport) model.

5.2.1.1 Groundwater Flow Field and Hydraulic Gradients

The modeller (or reviewer) should inspect the groundwater flow field to determine whether groundwater flow is predominantly in one, two, or three directions. Groundwater flow in many larger (regional) aquifer systems is predominantly horizontal with only a small vertical flow component. The relative contributions of horizontal versus vertical flow can be estimated using simple Darcy calculations (assuming appropriate estimates of horizontal/vertical hydraulic conductivity and hydraulic gradients). If the vertical flow component is small (say less than 5%) a two-dimensional representation of the regional groundwater flow field (using a plan view model) may be adequate.

Very few natural groundwater flow systems are one-dimensional due to the influence of local boundaries (e.g. valley boundaries, bedrock topography). Hence, a one-dimensional flow model is rarely appropriate, in particular for smaller aquifers which are strongly influenced by local boundaries. A one-dimensional approximation of groundwater flow may be adequate to simulate groundwater flow in a certain portion of a more complicated flow field. For example, groundwater flow from a relatively uniform hill side to a larger valley aquifer may be represented using a one-dimensional flow model.

In practice, the scale of groundwater resource projects and their influence on the groundwater flow field (see below) rarely justify the use of a one-dimensional flow model. One-dimensional models are predominantly used for initial scoping calculations (e.g. for conceptual modelling) typically in combination of simple analytical solutions such as a one-dimensional version of Darcy's Law (see below).

5.2.1.2 Degree of Heterogeneity

The modeller (or reviewer) should consider the degree of heterogeneity in the aquifer system when selecting the appropriate dimension of the mathematical model (see Section 4).

For example, the presence of a confining till layer separating a shallow aquifer potentially influenced by seepage from a mine waste unit from a deeper drinking water aquifer likely requires a discrete representation of this unit. In this scenario the potential movement of fluids between the surficial and deeper aquifers is a primary issue of concern. In order to model vertical flows, this scenario requires at least a two-dimensional cross-sectional model if flow in the transverse direction is small, or a full three-dimensional model.

Alternatively, highly permeable units which could represent preferential flow channels (e.g. a paleochannel or a fracture zone) may require explicit representation (in particular if contaminant transport is of concern) requiring a 2D or even a 3D modelling approach.

5.2.1.3 Type of Hydraulic Stresses and Influence on Flow Field

The modeller (or reviewer) should consider the type of stresses in the aquifer system when selecting the appropriate dimension of the mathematical model.

The most common stresses simulated in the natural resource industry include pumping, pit dewatering, underground dewatering and mine waste seepage (from TSF, WRDs and/or backfilled pit/underground workings.

Pumping from an extraction well creates a radial flow field which is best represented by a 3D axi-symmetric model or a 2D plan view model. A three-dimensional representation of a pumped aquifer may be required only if the pumping well penetrates a small portion of the aquifer ("partially penetrating" well) and/or the aquifer has significant vertical heterogeneity (e.g. confining layers).

Excavation of an open pit and associated dewatering tend to create a significant drawdown in the surrounding aquifer(s) somewhat analogous to a pumping well. In most mining projects, open pits reach significant depths and a representation of the vertical flow field is important. If the pit geometry is regular and the surrounding groundwater flow field is relatively uniform, a cross-sectional model may be adequate to simulate flow to the pit. In a more complex setting, a fully three-dimensional representation of the open pit and the surrounding aquifer may be required.

Figure 5-2 shows an example of a proposed open pit located in a complex setting where a three-dimensional flow model is appropriate. First, the open pit reaches a significant depth (300m) and intercepts bedrock units with greatly varying hydraulic conductivities. Second, the groundwater flow field is not uniform due to the presence of a lake on one side of the open pit and a steep mountain on the other side of the open pit. Finally, the open pit is intersected by several faults which are inferred to have a higher hydraulic conductivity than the local bedrock. The complexity of this setting prevents the use of a simplified cross-sectional model to simulate pit dewatering.

Figure 5-2: Example of a three-dimensional pit inflow problem. The upper panel shows a 2D plan view of the site, illustrating the proximity of the open pit to a nearby lake to the west and high relief to the east. The open pit is about 300m deep and intersects bed.

Mine waste units (such as TSFs or WRDs) tend to have a significant foot-print area and therefore commonly have a large effect on the groundwater flow field. As a result, the simulation of seepage from these mine waste units requires at least a 2D plan view model and often a three-dimensional model. The large size and regular shape of many mine waste units may result in a relatively uniform flow field (in x-y) which may be simulated using a cross-sectional model (see Appendix C1 for example Case Study 1).

5.2.1.4 Type of Remedial Work and Influence on Flow Field

The modeller (or reviewer) should consider the type of remedial work and its influence on the groundwater flow field in selecting the appropriate dimension of the mathematical model.

For example, an assessment of the efficacy of a shallow seepage collection drain and/or a partially penetrating cutoff wall requires at least a cross-sectional model or even a fully three-dimensional model. In contrast, an assessment of a fully penetrating drain and/or cutoff wall could be simulated in a two-dimensional (plan view model).

5.2.1.5 Required Model Accuracy & Data Uncertainty

The required model accuracy and data uncertainty influences the selection of an appropriate dimension for the flow and/or transport problem. In many applications the required model accuracy is not high enough to justify the use of a fully three-dimensional model. High model accuracy is often not warranted because of the lack of field data.

For example, if hydraulic conductivity estimates for the main aquifer influencing pit inflow have an uncertainty of two orders of magnitude, a 10% or 20% error introduced by ignoring the radial component of groundwater flow to an open pit may be acceptable. Often, analytical solutions can be used to estimate the error introduced by reducing the dimensionality of the flow problem.

Any simplification of a three-dimensional flow problem to a 2D (or even 1D) problem should be justified. An estimation of the resulting error should be provided and the potential implications for model predictions should be discussed.

5.2.1.6 Scale Issues

Most groundwater models recently developed to assess environmental impacts for larger mining projects used a basin-wide approach, i.e. a single three-dimensional model that encompasses the local watershed(s) is used to assess various aspects of the proposed project (see Table 1-1 in section 1).

These basin-wide groundwater models have the advantage that the domain boundaries are typically well-defined and that all three flow components are included. These 3D models are often not well-suited to evaluate specific aspects of groundwater flow due to the large scale of the model; however a coarse-grid regional model may be used to define boundary conditions for a fine-grid local site model.

As another example, the design of a passive pit dewatering system (with drain holes at various pit benches) may be better evaluated using a detailed cross-sectional model (using Cartesian or radial coordinates). This way, additional detail in the vertical (e.g. presence of pit benches, and confining layers) can be included in the 2D model that would unnecessarily "burden" the 3D basin model.

These guidelines recommend the use of a combination of modelling approaches from basin-wide models (in 3D) to detailed sub-models (ranging from 1D to 3D) to assess different aspects of the groundwater flow system.

Case Study 1 provides a good example of the use of multiple models for an assessment of environmental impacts for a proposed open pit mine (see Section 5.6.2).

5.2.2 Time Dependency

Almost all groundwater systems are inherently transient (dynamic) and the modeller has to determine whether this time dependency should be explicitly modeled or whether a steady-state (essentially does not change over time) simulation is adequate to satisfy the modelling objectives. Note that solute transport is, by definition, a transient process. Solute transport may be simulated assuming steady-state or transient flow conditions (see Section 9).

Many groundwater flow systems in British Columbia are highly transient, including:

• Aquifer(s) with strongly seasonal recharge (e.g. snowmelt dominated recharge in high elevation and/or cold region sites)
• Aquifer(s) influenced by variable pumping stresses (e.g. seasonal pumping, multiple users)
• Aquifer(s) influenced by seasonal irrigation (e.g. farmland)
• Coastal aquifer(s) or aquifers near large rivers which are tidally influenced (e.g. Fraser River).

The use of a transient groundwater flow model should be considered for any of these highly transient groundwater flow systems listed above.

A review of recent modelling studies carried out for EA submissions in British Columbia indicates that most flow models used for proposed mining projects assumed steady-state conditions (see Table 1-1 in section 1). Transient flow models were more commonly used in the assessment of groundwater resource projects which often require transient calibration of the model to pumping tests (Table 1-1).

The modeller has three options to represent a transient groundwater flow system:

• Simulate transient flow conditions (transient model)
• Simulate average flow conditions (steady-state model)
• Simulate high and low flow conditions (steady-state model).

Given the highly transient nature of most groundwater systems in B.C., and the complexity of a transient model, the first attempt to characterize the groundwater system should include a range of simulated conditions to bracket the existing variability, and the expected changes resulting from the proposed development.

It is good modelling practice to first develop and calibrate a steady-state model using average flow conditions for model calibration. If required, the transient aspects of the flow (or transport) problem can then be modeled using the calibrated steady-state model as initial conditions.

Fully transient models are an order of magnitude more difficult to develop and calibrate than steady-state models so careful consideration should be given to whether a transient model is required. Which modelling approach is selected depends on the flow system studied and the modelling objectives.

Consider the example of a flooded open pit which is located in an area with highly seasonal recharge (refer to Figure 4-3a/b in section 4). The flooded open pit is hydraulically connected to a bedrock aquifer which shows significant seasonal variations in groundwater levels (Figure 5-3). During the high-flow period, the groundwater is recharging the open pit and during the winter low-flow period, the open pit is recharging the aquifer. A steady-state model (using average recharge conditions) simulates the "average" water levels and the "average" (or net) flux to/from the open pit (Figure 5-3). However, a transient model (using the seasonal recharge distribution) would be required to assess the seasonal inflows and outflows to/from the open pit. In this example, the seasonal variations in (contaminated) groundwater discharge to the flooded open pit (and the nearby river) were important, and hence a transient model was required.

When simulating a transient problem assuming steady-state conditions, the modeller (or reviewer) should ensure that this approximation yields conservative model predictions. In the above example, high flow (steady-state) conditions could have been simulated to provide a conservative estimate of discharge of contaminated groundwater to the open pit and nearby rivers.

Any simplification of a transient flow problem to a steady-state problem should be justified. It should be demonstrated that this simplification yields conservative estimates of environmental impacts.

Figure 5-3: Comparison of steady-state and transient flow solutions for a highly transient groundwater flow system. Upper left panel shows estimated monthly recharge to the site and mean annual average. Lower left and upper right panels compares observed groundwater levels with simulated hydraulic heads (using transient and steady-state model). Lower right panel shows simulated groundwater flow from/to open pit using transient and steady-state model.

5.2.3 Type of Flow Domain

Three different approaches are available to simulate groundwater flow and solute transport in natural aquifer systems:

• Equivalent porous medium ("EPM")
• Discrete fracture network ("DFN")
• Dual porosity medium ("DPM").

The EPM approach assumes that the aquifer system can be represented by an equivalent porous medium, i.e. that the aquifer system behaves like a porous medium and standard flow and transport equations apply. The EPM approach is commonly used for unconsolidated materials such as overburden soils (colluvium), fluvial, alluvial and glacio-fluvial sediments, and highly weathered bedrock with high primary porosity.

The EPM approach is also commonly used to describe groundwater flow through fractured bedrock in which the primary porosity is very low and the effective permeability is controlled by fractures, fissures and bedding planes (i.e. secondary permeability). This approach is based on the assumption that at a sufficiently large scale (i.e. the representative elementary volume or REV) the bedrock mass will behave like a porous medium and can be described by "effective" hydraulic properties (e.g. specific discharge or flow per unit cross-sectional area).

The majority of groundwater modelling codes (including the codes discussed in Section 5.4.2) use the EPM approach to model groundwater flow.

In the discrete fracture network approach, it is assumed that flow through the bedrock matrix is negligible and all groundwater flow occurs through an interconnected network of fractures. Such a discrete fracture network may either be described explicitly (with known geometry) or generated randomly using fracture network statistics (e.g. Dershowitz et al., 2004; Parker and Cherry, 2011). Sophisticated modelling codes are available to generate DFNs and to simulate groundwater flow and solute transport in such a medium, including FracMan (available from http://www.fracman.com/ ) and Fractran (available from http://www.waterloohydrogeologic.com/software/fractran/fractran_ov.htm ).

Flow and transport in fractured bedrock and structured porous media (e.g. fractured sandstone) can be described using dual porosity models (DPM). This approach assumes that the medium consists of two regions, one associated with the macropore or fracture network and the other with a less permeable pore system of soil aggregates or rock matrix blocks (Gerke and van Genuchten, 1993). Different models exist to describe the nature of flow and transport in these two domains and the extent of their interaction. In its simplest form, groundwater flow and advective transport is assumed to only occur in the highly permeable ("active") domain. Groundwater flow in the low-permeable ("inactive") domain is assumed to be negligible but this stagnant zone influences solute transport by diffusion.

At present the DFN and dual porosity models are predominantly used in research and/or in assessment of contaminated sites with very high risk and/or consequence (e.g. storage of radionuclides, large contaminated sites impacting drinking water supplies, etc.). The primary challenge with the DFN and DPM models is model parameterization. A characterization of the fracture network and/or the dual porosity regime requires extensive field studies and/or detailed model calibration usually not available for natural resource projects.

These guidelines generally recommend the use of the EPM approach for groundwater modelling in support of natural resource projects in fractured bedrock settings but the application must be reasonably justified (e.g. supporting site data, simple objectives, etc.). Note, however, that the use of the EPM approach in moderately to sparsely fractured bedrock significantly increases model uncertainty. For example, an EPM model may correctly predict the overall inflow to an open pit but may not be able to predict the exact location where such an inflow will occur. Even greater uncertainty should be expected when modelling transport processes with an EPM approach, as the direction and timing of contaminant transport can be strongly influenced by flow through discrete fractures (see also section 9). If applied to fractured bedrock, the limitations of the EPM approach and potential implications for model predictions should be discussed in the model report.

DFN or DPM approaches should be considered in projects where the modelling objectives justify these approaches (e.g. contaminant transport in fractured bedrock with high risk/consequence) and the supporting site characterization data is available.

5.2.4 Variably Saturated Flow

Most traditional groundwater models simulate "saturated flow", i.e. water movement in the saturated zone below the water table. In this approach, water movement in the "vadose zone", i.e. the unsaturated zone above the water table, is ignored. For example, most groundwater models apply recharge directly to the top of the aquifer without simulating the processes of infiltration through the vadose zone, transpiration from the root zone, etc.

During model scoping, the modeller (or reviewer) should determine whether water movement in the vadose zone should be explicitly simulated. Water movement in the vadose zone is described by the Richards equation (Richards, 1931) and requires knowledge about the unsaturated hydraulic conductivity which is a function of the degree of saturation of the medium.

In the context of natural resource projects, an in-depth assessment of water movement in the vadose zone may be required for the following conditions:

• Assessment of seepage from a mine waste facility (e.g. a tailings impoundment)
• Simulation of groundwater flow with significant influence of perched zones
• In-depth assessment of soil cover performance.

Several "variably saturated" flow models have been developed which simulate the continuum from unsaturated flow in the vadose zone to saturated flow in the saturated aquifer (Section 5.4.2). Solution to the Richards equation is non-linear and requires significantly greater computing effort than conventional saturated flow models. This approach requires significantly greater vertical discretization to achieve model convergence.

For the reasons listed above, variably-saturated flow (and transport) models are commonly simplified to 2D or even 1D flow problems to reduce model complexity and computing time.

The use of variably saturated flow for large-scale models of groundwater flow (e.g. at the basin scale) is uncommon and not recommended. In most cases inadequate field data are available to either parameterize and/or calibrate the model. This tends to result in over-parameterization which makes the model less transparent and unnecessarily complex.

In some instances, a variably-saturated flow solution may be used to avoid numerical problems associated with simulation of the water table ("phreatic surface") in regional aquifers (Section 5.4.2). Some modelling codes provide options for simplified formulations of the hydraulic conductivity function to use this approach at the regional scale (e.g. FEFLOW; MODFLOW-SURFACT). If this approach is used, the model parameters often do not represent realistic physical parameters (because the simplified formulations are not physically based).

If a variably-saturated flow solution is used to solve a 3D flow problem, this approach and the hydraulic parameters assumed for the unsaturated zone, should be justified and documented.

5.2.5 Solute Transport Model versus Particle Tracking

If the modelling objectives for a given project include an assessment of contaminant transport, then the modeller (or reviewer) should determine whether solute transport can be simulated using particle tracking or whether a full solute transport model is required.

Particle tracking visualizes the flow path of a solute and allows an estimation of travel times assuming advective transport. This method does not allow the prediction of contaminant plume migration in space and time (see Section 9 for more details).

Solute transport models simulate the movement of contaminant plumes assuming all dominant transport processes (i.e. advection, dispersion, and diffusion). These models allow a quantitative assessment of contaminant transport, including the prediction of the spatial distribution of contaminant concentrations in the aquifer and/or contaminant breakthrough at a given location (see Section 9 for more details).

Particle tracking provides significant insight into the groundwater flow field and the resulting flow paths and should always be used, at least as a screening tool, for contaminant transport problems.

Solute transport modelling should be considered for an impact assessment of natural resource projects if potential CoCs have been identified in the project AND any of the following conditions apply:

• Particle tracking suggests that CoCs may reach VEC(s) (e.g. streams) via groundwater, and
• Conservative mixing calculations (i.e. mixing of contaminant load into receiving water without chemical reactions based on particle tracking and/or load balance modeling) suggest that contaminant concentrations may exceed applicable guidelines
• Significant changes in the source concentration of a CoC are expected over the course of the project
• Significant dispersion of the contaminant plume is anticipated and poses a risk to the environment
• Significant dilution of the contaminant of concern is expected along the flow path due to dispersion and/or dilution (e.g. by recharge or groundwater inflow); or
• The CoCs are known to be reactive and geochemical reactions may significantly influence contaminant transport (see section 5.2.6 below)

These guidelines recommend the use of a phased approach for assessing contaminant transport and associated water quality impacts for natural resource projects. Initially, particle tracking should be used to determine the potential flow paths and the risk to VECs. Preliminary (conservative) estimates of contaminant loads and/or contaminant concentrations to VECs (e.g. rivers, lakes) should be obtained to assess whether applicable guidelines may be exceeded. If these conservative calculations indicate a potential risk to VECs, then solute transport modelling should be considered.

Note that solute transport modelling may not be required if the project descriptions includes mitigation options (e.g. lining of a tailings facility, placement of waste rock under water) that eliminates or at least adequately reduces the contaminant source(s).

5.2.6 Conservative versus Reactive Transport

Solute transport in the subsurface is commonly influenced by chemical reactions occurring: (i) among different solutes (dissolved in groundwater); and/or (ii) with the surrounding solid phase (i.e. the soil particles or rock surface).

Common reactive transport processes encountered in aquifers include: (i) sorption/desorption; (ii) precipitation/dissolution; and (iii) redox reactions (see Section 9 for more detail). Another category of reactive processes is biologically facilitated transformations, like dentrification. Although reactive transport processes have been documented and studied for many years, their influence can vary greatly from site to site. As a result, site-specific information is usually required to quantify their effects on contaminant transport.

In most cases, reactive transport processes tend to reduce contaminant concentrations. Hence, consideration of reactive transport processes tends to yield non-conservative water quality predictions as there are sinks for the solute. For this reason, the use of reactive transport models should be limited to situations where all of the following situations apply:

• Conservative transport modelling has demonstrated that there is a potential exceedance of applicable standards in one or several CoCs.
• Reactive transport processes included in the model are well-established in the literature.
• Reactive transport parameters used in the model are
  • based on site-specific laboratory/field testing and/or
  • calibrated against field observations (e.g. breakthrough or plume distribution.)

The above guidelines are consistent with the phased approach for contaminant transport described recommended in Section 5.2.5. Conservative transport modelling should precede reactive transport modelling to identify whether the additional complexity of reactive transport modelling is justified. This phased approach ensures that the influence of the reactive transport assumptions on water quality predictions is transparent.

If a reactive transport model is selected, the modeller should justify this decision, and provide details on model parameterization and calibration, including supporting lab/field studies.

Note again that solute transport modelling (whether conservative or reactive) may not be required provided the project description includes mitigation options (e.g. lining of a tailings facility, placement of waste rock under water) that would eliminate (or adequately reduce) the risk to any VECs (and this can be demonstrated using conservative mass balance calculations and/or particle tracking.

5.3 Overview of Modelling Tools

The groundwater modeller has a wide variety of modelling tools available at his/her disposal to model the groundwater system. The modelling tools can be broadly grouped into three categories:

• Analytical Models
• Numerical Models
• Analytic Element Models.

This section provides an overview of these modelling methods and describes their main advantages and disadvantages.

5.3.1 Analytical Models

5.3.1.1 Definition

Analytical models use exact solutions to the equations that describe groundwater flow or contaminant transport. In order to produce these exact solutions, the flow/transport equations have to be considerably simplified such that they are typically applicable only to simple flow and contaminant transport systems. Analytical models can be simple formulae, spreadsheets, or sequences of calculations packaged in a piece of software.

5.3.1.2 Model Use

Table 5-1 lists the advantages and disadvantages of analytical model vis-à-vis numerical models. The main advantage of analytical models is the ease of use and transparency of such models which will facilitate sensitivity analyses. Their main disadvantage is that they can only be applied to relatively simple flow (or transport) problems.

The main uses of analytical models are to:

• Assist in conceptual modelling
• Simulate flow and/or transport in simple physical settings (or where there are only one or two simple objectives)
• Simulate flow and/or transport for projects which have low risk of impact to VECs via groundwater
• Check results of the numerical model.

Analytical models should be used as a starting point in the assessment (e.g. during conceptual modelling) and before moving on to more sophisticated numerical models. Experimentation with analytical equations to examine the influence of changing parameter values on the model results is a vital exercise to gain an understanding of which are the key (i.e. the most sensitive) parameters and to develop a correct conceptual model.

For projects with relatively simple physical settings, or where the required accuracy of model predictions is not very high, analytical models may be adequate to simulate groundwater flow and/or contaminant transport. Analytical models may be the preferred choice of modelling for projects where there is insufficient monitoring data available, i.e. where the uncertainty due to data limitation is greater than the error caused by using a simplified analytical model. However, a lack of data does not justify reliance on an analytical solution alone; for complex systems, more data is required to support a more thorough modelling effort.

Where there are only one or two simple modelling objectives (i.e. determining drawdown or flow to a well), an analytical model may be adopted; however, where more modelling objectives are identified, a numerical model may provide more versatility.

Table 5-1: Advantages and Disadvantages of Analytical and Numerical Methods.

The use of analytical solutions may also be adequate for projects with low risk of impact to VECs via groundwater. This includes projects where the project description includes mitigation options that reduce or even eliminate the contaminant source to such an extent that there is no (or only a very low) risk to VECs.

Analytical models may also be used by the modeller (or the reviewer) to check that numerical solutions produce approximately the right answer. It is very easy to mistype input for complex models and an analytic check can generate confidence in the numbers produced by a numerical model.

Note that analytical models may be used in conjunction with numerical models. For example, analytical models can be used to run preliminary sensitivity analyses on the groundwater system to define the most appropriate parameters for inclusion in a steady-state model. Furthermore, analytical transport solutions (e.g. Ogata Banks equation) may be used in conjunction with a numerical groundwater flow model. In this situation, the analytical transport equation could be used to estimate the breakthrough of a contaminant along an identified flow path which was determined by the numerical flow solution. The use of an analytical transport solution greatly reduces the computing time, in particular, if a detailed sensitivity analysis of various transport parameters is required.

5.3.1.3 Application to Natural Resource Projects

Analytical models are available for a wide range of flow and transport problems which are applicable to the natural resource industry, including:

• Groundwater flow in 1D and 2D (e.g. for conceptual modelling)
• Groundwater flow to a well (e.g. for pump test analysis, pit inflow estimates.)
• Groundwater flow to a trench (e.g. for inflow to collection ditches, underground mine.)
• Groundwater flow to a tunnel (e.g. for inflow to adits, underground workings.)
• Groundwater flow to an open pit (e.g. for inflow to open pit mine)
• Solute transport with 1D flow (e.g. for transport modelling along flow path).

Table 5-2 lists some useful analytical solutions commonly used in the resource industry. For more details on analytical solutions the reader is referred to Bear (1972), Domenico and Schwartz (1990) and Fetter (1992). A compilation of analytical solutions to common transport problems is provided in a software program (STANMOD) which was developed by the US Salinity Laboratory (Simunek et al., 1999) and can be downloaded from their website at http://www.ars.usda.gov/Services/ .

Table 5-2: Useful analytical solutions for the resource industry.

5.3.2 Numerical Models

5.3.2.1 Definition

A numerical model uses numerical methods to solve the governing equations of groundwater flow and/or contaminant transport. In distributed numerical models, space and time are divided into discrete intervals (as illustrated by Figure 5-4 and 5-5) where for each model grid cell, parameter values are defined including hydraulic conductivity, porosity, aquifer thickness, initial contaminant concentration, etc.

Numerical models enable more complex systems to be represented than can be represented by analytical models. Furthermore, numerical models may allow for multiple modelling objectives to be addressed in parallel. Numerical models still require simplifications to be made about system behaviour.

Figure 5-4: Illustration of finite-difference approach (reproduced from NGCLC, 2001).

Figure 5-5: Illustration of finite-element approach (reproduced from NGCLC, 2001).

5.3.2.2 Model Use

Table 5-1 lists the advantages and disadvantages of analytical models vis-à-vis numerical models. The main advantage of numerical models is that different parameter values can be assigned to each cell, so that lateral and vertical variations in property values can be taken into account. The geometry of the model can be designed to reflect the geometry of the system. In addition, models can be constructed that include more than one layer; this enables multi-layered aquifers to be represented. For time variant models, model inflows (e.g. recharge and its contaminant concentration) and outflows (e.g. groundwater abstractions) can be specified for each model time step.

Their main disadvantage is that numerical models can be costly and time-consuming. Another potential disadvantage is that the model complexity reduces the transparency of the model calculations and/or can mask the remaining model uncertainty.

Numerical models will generally be applicable where:

• Previous modelling studies (or conceptual modelling) using simple analytical models have shown that a more sophisticated approach, such as incorporating spatial variability, is required.
• The groundwater regime is too complex to be robustly represented by an analytical model.
• The required model accuracy (as defined by the model objectives) requires the use of a numerical model.
• Processes affecting contaminant transport cannot be adequately represented by simple transport equations.
• An analytical model is inadequate for the design of mitigation measures, e.g. in determining the optimal location and pumping rate for boreholes in a pump and treat scheme.

Numerical models should be considered where the scale and importance of the problem warrant the use of a more sophisticated approach. For such sites, the scale of the problem should demand detailed site investigations which should provide sufficient information to allow the construction of a numerical model.

Numerical models should not be used as an alternative to data collection. Instead, appropriate use of a numerical model may require and/or help guide additional data collection so that the model can be properly parameterized and calibrated as well as support system interpretation.

5.3.2.3 Application to Natural Resource Projects

Numerical models are versatile and flexible with respect to model domain, boundary conditions and other stresses and can be used to simulate almost any groundwater flow and contaminant transport problem related to natural resource projects in British Columbia (See section 3).

A range of different numerical models (also referred to as "modelling codes" or simply "codes") are available to cover specific aspects of groundwater flow and contaminant transport (see Section 5.4 below).

Some numerical models only solve the equations for groundwater flow ("flow models") while other models solve only the equations for contaminant transport ("transport models") and yet others solve both sets of equations ("flow and transport models"). Some modelling codes are very versatile while other modelling codes have a narrower range of application (e.g. 2D cross-sectional models) or specialize on certain aspects of groundwater flow (e.g. variably saturated flow, density-dependent flow).

A more detailed review of modelling codes suitable for use in the resource industry is given in Section 5.4 below.

5.3.3 Analytic Element Models

5.3.3.1 Definition

An analytic element model uses superposition of closed-form (analytical) solutions to the governing differential equation of groundwater flow to approximate both local and (near-field) and regional (far-field) flow. Hence, analytic element models do not require grid discretization or specifications of boundary conditions on the grid perimeter (Hunt et al., 1998).

These characteristics allow for representation of large domains that include many hydrogeologic features outside the immediate area of interest (i.e. far-field) and easy modification of the regional flow field by adding analytic elements representing regional hydrologic features.

5.3.3.2 Model Use

Analytic element models are well-suited for use as screening models (Hunt et al., 1998). Analytic element models can be used to develop conditions on the grid perimeter for a smaller numerical model, similar to the process of telescopic mesh refinement (TMR). The advantage over traditional TMR using finite difference models is that this method: (i) allows easy addition of far-field elements until the far field is correctly simulated; and (ii) avoids discretization problems that can occur in large-scale models with large cell/element sizes.

The major limitation of analytic element models is that the method is only computationally efficient for steady-state flow in large aquifers where the vertical flow component can be ignored.

5.3.3.3 Application to Natural Resource Projects

The analytic element method is not commonly used by the consulting community to study natural resource projects, due to lack of easy-to-use software codes, familiarity by the model practitioners, and/or its limitations (only applicable to 2D steady-state flow).

5.4 Numerical Modelling Tools

5.4.1 Solution Methods

Several numerical methods are available to solve the governing equations for groundwater flow and contaminant transport, including (NGCLC, 2001):

• Finite difference method (flow & transport)
• Finite element method (flow & transport)
• Eulerian methods such as TVD (transport)
• Lagrangian method (transport)
• Mixed Eulerian-Lagrangian method (transport).

The following sections provide brief overview of the two common methods to solve groundwater flow and transport. The solution methods available to specifically solve solute transport are discussed further in Section 9.

5.4.1.1 Finite Difference Method

Finite Difference Method ("FDM") is the most commonly used approach in numerical groundwater flow modelling. For most finite difference models the space and time co-ordinates are divided on a rectangular grid (Figure 5-4), and model parameters (such as hydraulic conductivity, aquifer thickness etc.) are specified for each model grid cell. The flow and transport equations are solved by direct approximation. The grid spacing represents the degree of accuracy of the model in representing lateral or vertical changes in the property values that describe the system. Finite difference methods have the advantage of being relatively simple to use, but have the disadvantage of not accurately representing irregular boundaries; also it is difficult to change the grid spacing to provide greater precision in areas of interest.

5.4.1.2 Finite Element Method

In the Finite Element Method ("FEM") the spatial domain is divided into a mesh of elements, generally of triangular or quadrilateral shape (Figure 5-5). Variation in a model parameter across the model element is normally approximated by a polynomial function. This technique provides greater flexibility than finite difference methods in representing the model domain, particularly complex geological boundaries. The model mesh can easily be modified to provide greater precision in areas of interest, although complex meshes require software tools for their management. Finite element models are less susceptible to numerical dispersion than finite difference models, but for the same number of elements/cells the computing requirements are higher.

5.4.2 Model Codes

A groundwater modelling code is defined as a computer code that solves a groundwater flow or contaminant transport problem. The computer code facilitates input of relevant model input parameters (e.g., model grid, aquifer parameters, boundary conditions.) and solves the groundwater flow or contaminant transport problem. Most model codes provide options for viewing the model output (e.g. head solutions in space and time; groundwater budget etc.).

Table 5-3 lists commercially available groundwater modelling codes commonly used in the groundwater industry in North America. The following codes are routinely used in the resource industry in British Columbia and are discussed below in more detail:

• MODFLOW (flow)
• MODFLOW SURFACT (flow and transport)
• MODPATH (particle tracking)
• MT3DMS (transport)
• SEEP/W (flow)
• FEFLOW (flow & transport).

The list of modelling codes described in Table 5-3 and discussed here does not include all available modelling codes. Many other groundwater modelling codes exist and new groundwater modelling codes will undoubtedly be developed in the future. The groundwater modeller is encouraged to stay abreast of the development of new computer codes to find the code most suitable to meet the project objectives.

Table 5-3: Commonly used numerical groundwater modelling codes in North America (adapted from MDBC, 2001).

5.4.2.1 MODFLOW

MODFLOW is a three-dimensional finite-difference ground-water model that was developed by the USGS and first published in 1984 (McDonald and Harbaugh, 1984). This code has a modular structure that allows it to be easily modified to adapt the code for a particular application. Many new capabilities have been added to the original model. The most recent version is MODFLOW-2005 (Harbaugh, 2005) but the earlier version MODFLOW-2000 (Harbaugh et al., 2000) is still in common use.

MODFLOW simulates steady and non-steady ("transient") flow in an irregularly shaped flow system in which aquifer layers can be confined, unconfined, or a combination of confined and unconfined ("convertible layer"). Flow from external stresses, such as flow to wells, areal recharge, evapotranspiration, flow to drains, and flow through river beds, can be simulated. Hydraulic conductivities or transmissivities for any layer may differ spatially and be anisotropic (restricted to having the principal directions aligned with the grid axes), and the storage coefficient may be heterogeneous. Specified head and specified flux boundaries can be simulated as can a head dependent flux across the model's outer boundary.

MODFLOW has a built-in solute transport function (MOC-3D). MODFLOW is more commonly used in combination with MT3D, another modular solute transport code.

MODFLOW is the most commonly used groundwater code world-wide and many graphical user-interfaces (GUI) have been developed to support its use (see below). This code is well-documented and verified.

The main advantages of MODFLOW are that the code is robust, easy to use, and versatile. In addition, the availability of powerful GUIs greatly facilitates the model setup (pre-processing) as well as interpretation and visualization of model output (post-processing). The MODFLOW code also enables the modeller to extract detailed water balance information from the model (using the FLOWBUDGET subroutine) which greatly assists with model interpretation and trouble shooting.

The main limitations of MODFLOW include: (i) the inefficient use of rectangular grids to represent complex geometries; and (ii) convergence problems related to wetting/drying of cells near the water table. The first limitation is common to all finite-difference models and can be dealt with by using a denser grid spacing (which is no longer a major limitation with increased processing power).

The second limitation is more important, in particular for unconfined flow problems in steep terrain. When solving an unconfined flow problem, MODFLOW checks the saturated thickness of the uppermost (active) unconfined (or convertible) layer. If the simulated water table in a given cell falls below the bottom elevation of the cell, this cell is made inactive (or "dry"), In subsequent iterations, this cell may be made "active" again (i.e. "wetted"). This wetting capability can cause significant convergence problems, in particular in steep terrain where the water table may interact with several model layers.

Note that the USGS has recently released a new MODFLOW code that can handle the drying/wetting issue (MODFLOW-NWT). At present, this version is still not included in most GUIs for MODFLOW and is therefore not widely used.

In many MODFLOW applications to natural resource projects (where steep topography is common) this problem is avoided by assuming confined conditions. This simplification may be adequate for some situations (e.g. baseline conditions) but can result in significant error in flow predictions (e.g. in inflow to an open pit).

These guidelines generally recommend the use of MODFLOW for most natural resource projects. More sophisticated codes (such as MODFLOW SURFACT or FEFLOW) may be required for projects with complex geometry and/or in steep topography (see below).

5.4.2.2 MODFLOW SURFACT

MODFLOW-SURFACT is a proprietary code developed by Hydrogeologic Inc. to simulate saturated/unsaturated groundwater flow and solute transport (HGL, 2008). MODFLOW-SURFACT was developed to overcome numerical difficulties encountered with MODFLOW, primarily related to the drying/wetting problem.

MODFLOW-SURFACT solves the fully 3D saturated/unsaturated groundwater flow equations, or alternatively, solves enhanced equations for performing unconfined simulations to rigorously model desaturation/resaturation of aquifers.

Additional improvements offered by the MODFLOW-SURFACT code include:

• Use of a curvilinear grid for efficiently fitting irregular domain geometries
• Additional boundary conditions (seepage face, unconfined recharge)
• Adaptive time stepping and restart options.

The primary advantage of MODFLOW-SURFACT is the handling of complete desaturation and resaturation of grid blocks and accurate delineation and tracking of the water table position. These additional abilities may be of importance in some natural resource projects, in particular where unconfined conditions are encountered in steep terrain (or is induced by mine dewatering).

These guidelines recommend the use of MODFLOW-SURFACT for natural resource projects where the vertical component of unconfined groundwater flow is an important modelling objective.

5.4.2.3 MODPATH

MODPATH is a particle-tracking post-processing package that was developed to compute three-dimensional flow paths using output from steady-state or transient ground-water flow simulations by MODFLOW. MODPATH uses a semi-analytical particle tracking scheme that allows an analytical expression of the particle's flow path to be obtained within each finite-difference grid cell.

MODPATH calculates pathlines and travel times of groundwater flow (and solutes dissolved in groundwater). Particle tracking can be used to draw flow nets, determine recharge zones and/or capture zones, and to estimate travel times of conservative contaminants of concern (see Section 9).

MODPATH is the most common particle tracking code in the groundwater industry. It is useful as a visualization tool to help understand flow patterns in simulated ground-water flow systems. It is useful for delineating sources of water to discharge sites and aquifers in systems simulated with MODFLOW.

This code is efficient, very easy to use and provides good visualization options (e.g. time markers). These guidelines recommend the use of MODPATH in conjunction with MODFLOW for particle tracking.

5.4.2.4 MT3DMS

MT3DMS is a comprehensive three-dimensional numerical model for simulating solute transport in complex hydrogeologic settings (Zheng and Wang, 1999). MT3DMS has a modular design that permits simulation of transport processes independently or jointly. MT3DMS is capable of modelling advection in complex steady-state and transient flow fields, anisotropic dispersion, first-order decay and production reactions, and linear and nonlinear sorption.

MT3DMS is linked with the USGS groundwater flow simulator MODFLOW, and is designed specifically to handle advectively-dominated transport problems without the need to construct refined models specifically for solute transport.

MT3DMS combines three major classes of transport solution techniques in a single code, namely, the standard finite difference method; the particle tracking based Eulerian-Lagrangian methods; and the higher-order finite-volume total-variation-diminishing (TVD) method. This unique range of solution techniques allows the user to solve a wide variety of transport problems ranging from advection-dominated (e.g. using TVD) to mixed advection-dispersion problems (EL methods) to dispersion dominated problems (FD).

MT3DMS is implemented with an optional dual-domain formulation for modelling mass transport in highly heterogeneous porous media or fractured media with a mobile domain (where solutes are moved by groundwater flow) and an immobile domain (where no groundwater flow occurs and solutes only move by diffusion).

This code is very well documented and is supported by all major GUIs developed for MODFLOW.

These guidelines recommend the use of MT3DMS for the simulation of solute transport problems of natural resource projects (in conjunction with MODFLOW).

5.4.2.5 SEEP/W and VADOSE/W

SEEP/W is a finite-element code for analysis of 2D seepage and excess pore-water pressure dissipation problems in porous media (Geoslope, 2007). SEEP/W can simulate steady-state confined and unconfined flow, transient flow, 2-D flow in a cross-section or in plan view, and 3D axisymmetric flow.

SEEP/W can simulate both saturated and unsaturated flow, a feature that greatly broadens the range of problems that can be analyzed. In addition to traditional steady-state saturated flow analysis, the saturated/unsaturated formulation of SEEP/W makes it possible to analyze seepage as a function of time and to consider such processes as the infiltration of precipitation. The transient feature allows you to analyze such problems as the migration of a wetting front and the dissipation of excess pore-water pressure.

Boundary condition types available in SEEP/W include total head, pressure head, or flux specified as a constant or a function of time; transient flux as a function of computed head; and review and adjustment of seepage face conditions.

SEEP/W is commonly used for engineering problems, in particular for assessment of seepage and associated pore pressures in tailings dams. However, this code is also capable of simulating variably saturated flow in natural aquifer systems (e.g. perched layers), but only in a 2D cross-sectional domain. Therefore, seepage from the lens in the transverse direction would need to be negligible, which may not be reasonable for smaller lenses.

SEEP/W has an easy-to-use graphical user interface (GUI) which facilitates setup of the model (mesh construction, model parameterization) and visualization of the model output (flow net, pore pressures etc.).

VADOSE/W has similar capabilities to those of SEEP/W but has a more rigorous algorithm to describe the soil-atmosphere interactions at the ground surface, including the simulation of actual evapotranspiration and net infiltration using atmospheric inputs.

These guidelines recommend the use of SEEP/W and VADOSE/W for the assessment of 2D seepage problems involving engineering structures (e.g. tailings dams) and 2D infiltration problems, respectively.

5.4.2.6 FEFLOW

FEFLOW is a versatile finite-element code that solves the following groundwater flow and contaminant transport problems (DHI-WASY, 2007):

• Transient or steady-state flow (3D)
• Saturated and unsaturated flow
• Multiple free surfaces (perched water table)
• Density-dependent flow (salt water intrusion)
• Contaminant and heat transport.

FEFLOW has a comprehensive selection of graphical tools for building the finite element mesh, assigning property zones and setting boundary conditions. The modelling platform includes state-of-the-art 3-D visualization of modelling outputs.

The main advantages of FEFLOW include: (i) versatility of code to solve different flow and transport problems; (ii) flexibility in model discretization due to use of the finite element method; and (iii) flexibility in formulation of boundary conditions. This flexibility in model discretization and boundary conditions can be very useful when simulating complex flow problems (e.g. mine developments in structurally controlled bedrock; progressive excavation of an open pit and/or underground mine workings).

FEFLOW can solve groundwater flow in the following modes:

• Unconfined ("Phreatic") mode
• Confined mode
• Variably saturated mode.

The use of the phreatic or the variably saturated option can lead to numerical instability and/or non-convergence if the water table crosses model layers (a common problem at project sites with steep topography). The assumption of a confined aquifer may be required to get a stable solution for problems with steep terrain (similar to MODFLOW).

Other limitations of the FEFLOW code which the modeller (or reviewer) should be aware of include:

• FEFLOW has limited capabilities to evaluate the water balance (e.g. change in storage is not provided in the water balance; flux section tool provides only approximate internal fluxes).
• FEFLOW does not allow pinching out of model layers and will simulate "artificial" flow through layers above the water table in phreatic mode (using the saturated hydraulic conductivity).
• Phreatic conditions may produce artificial water if residual water depth in "dry elements" is set too large.
• The transport algorithm of FEFLOW is prone to numerical oscillations and/or numerical dispersion requiring a high degree of horizontal and/or vertical discretization (which can be prohibitive for regional models).

The use of FEFLOW requires significant modelling expertise, including an in-depth understanding of finite-element methods and the FEFLOW code itself. The flexibility of the FEFLOW code makes it a powerful tool but a difficult one to use. For this reason, the use of FEFLOW is only recommended for experienced modellers.

These guidelines recommend the use of FEFLOW for complex natural resource problems in which complex geometries and/or complex boundary conditions will have to be simulated. FEFLOW is suitable for density-dependent groundwater flow problems (e.g. saltwater intrusion problems and/or deep groundwater flow involving brines).

5.4.2.7 Other Codes

As mentioned at the beginning of this section, the groundwater modelling codes listed in Table 5-3 and described above are those most commonly used by the industry in British Columbia. This list is not meant to be prescriptive and other groundwater modelling codes are available and may be equally, or better, suited for a given groundwater flow and solute transport problem.

If other, less common modelling codes are used, the user should confirm and document that the code in question has been adequately verified (see Section 5.4.4 below).

5.4.3 Graphical User Interface (GUI)

Most proprietary modelling codes come with a built-in Graphical User Interface (GUI) that improves model parameterization (pre-processing) and visualization of model output (post-processing).

Some modelling codes, namely MODFLOW and MT3DMS, include commercial graphical user interfaces that are commercially available. Table 5-4 summarizes the most common GUIs for MODFLOW/MT3DMS.

The capabilities and limitations of a GUI can significantly influence the outcome of a modelling study and should therefore also be considered when selecting a mathematical code for a given modelling objective. For example, some GUIs allow the use of automated parameter estimation while others do not support this option, thus limiting the options for model calibration.

Table 5-4: Common Graphical User Interfaces (GUI) for MODFLOW/MT3DMS (adapted from MDBC, 2001).

The following aspects of a GUI should be considered by the modeller during model selection:

• General
  • Reliability and stability of GUI
  • Technical support (User support and product development)
  • Flexibility in solution techniques (especially important for solute transport)
  • Support of automated parameter estimation (e.g. PEST).
• Pre-processing capabilities
  • Compatibility with GIS and/or ACAD
  • Ease of definition (and later adjustment) of boundary conditions
  • Ease of model discretization (and later regridding/remeshing)
  • Ease of implementation of complex geometries (e.g. boundary conditions along arcs,
• Post-processing capabilities
  • Support for model calibration (water balance output, visualization of residual error in scatter plots, histograms, maps etc.)
  • Ease & automation of modelling outputs (e.g. to Excel).

The selection of the most appropriate GUI depends on the project specifics and the modelling objectives. The capabilities and limitations of the GUI should be carefully considered by the modeller at the outset of the project. If limitations in the GUI are encountered at a later stage during the modelling study, the modeller should be prepared to switch to another GUI (and/or mathematical code).

These guidelines recommend the use of flexible GUIs which enable the user to easily adjust the conceptual model (domain boundaries, boundary conditions, zonation) and regrid/remesh the domain (e.g. GMS). This flexibility helps the modeller to start with a simpler model and gradually build-up the complexity of the model (rather than starting out with a very complex model).

5.4.4 Code Verification

Groundwater modelling codes used for the assessment of potential impacts by natural resource projects should be verified.

The process of code verification involves a check that the code is free of errors, i.e. mathematical equations are appropriately coded and the modelling results are correctly output. At a minimum, code verification should include a comparison of the modelling results against analytical solutions (for simplified problems). In addition, model verification may include a comparison of the modelling results against a numerical solution obtained using other (preferably verified) modelling codes (for more difficult problems).

All modelling codes (and supporting GUI) listed in Table 5-3 and described above have been extensively verified. However, continuous upgrading of modelling codes and/or improvements of the GUIs can introduce errors that may not be readily apparent. If numerical problems are suspected or encountered then selected modelling results (for example the calibrated model) could be checked using different modelling codes (e.g. MODFLOW and FEFLOW). A good agreement of modelling results using independent modelling codes will ensure that modelling results are not influenced by input errors and/or code selection.

The modeller should document all code verification completed as part of the QA/QC of the modelling study.

5.5 Code Selection Process

Figure 5-6 summarizes the process of selecting the appropriate mathematical model for groundwater flow and solute transport. The first step in the process is to determine whether the groundwater flow regime can be represented by an analytical model or whether a numerical model is required.

If an analytical solution is not available (or not sufficiently accurate) a numerical flow code should be selected. During the selection the capabilities and limitations of the available flow codes (see Section 5.4) should be compared to the defined modelling scope (see Section 5.2) and the code most suitable should be selected. In many instances, alternative codes are available to solve the same flow problem. In this case the final code selection may be influenced by other factors such as the availability of the code and/or the familiarity of the modeller with the code and GUI.

A similar selection process is then repeated for the contaminant transport problem (Figure 5-6). The first step in this selection process is to determine whether an analytical solution exists for the transport problem. Note that the use of a numerical flow model does not necessarily require the use of a numerical transport model. For example, a complex three-dimensional flow model may provide for the definition of a (simplified) 1D flow path. An analytical transport solution may then be applied to simulate transport along this 1D flow path assuming average hydraulic properties.

If no analytical transport solution is available (or not sufficiently accurate), a numerical transport code should be selected. Again, the capabilities and limitations of the available transport codes (see Section 5.4) should be compared to the defined modelling scope (see section 5.2) and the code most suitable should be selected.

As discussed in Sections 5.2.5 and 5.2.6, these guidelines recommend a phased approach to transport modelling which may require the use of different transport models, ranging from simple (conservative) particle tracking codes to sophisticated solute transport modelling codes.

Figure 5-6: Flow chart illustrating selection of mathematical model.

5.6 Case Studies

5.6.1 Case Study 1: Open Pit Mine

5.6.1.1 Overview

Case Study 1 illustrates the use of multiple mathematical models (and codes) to assess the potential impacts of a proposed open pit mine. A more detailed description of this case study is provided in Appendix C.

The following models were developed to address modelling objectives:

• A 3D basin-wide groundwater flow model to assess regional groundwater flow and pit inflow
• An analytical model (Darcy calculations) to estimate pit inflow along a fault
• A 2D (cross-sectional) variably-saturated flow model to assess seepage from the TSF
• A 3D local flow model to evaluate seepage from the TSF to local creeks (VEC).

The use of multiple models to address various groundwater flow aspects of the project (which differ in scale and complexity) is good modelling practice and is generally recommended.

5.6.1.2 Basin-wide model

Regional groundwater flow was conceptualized to be three-dimensional. The model domain covers an area of about 40 km² and is bounded primarily by the natural groundwater divides formed by creeks to the north, south, and east, and by bedrock topography to the west (Figure 5-7). Vertically, the model was discretized into five hydrostratigraphic layers representing glacial outwash, till, and bedrock.

The three-dimensional finite difference code MODFLOW-SURFACT was selected for the basin-wide modelling for two primary reasons:

• The code integrates groundwater and surface water flow systems, and;
• The code employs an updated wetting/drying function which minimizes convergence problems typically associated with MODFLOW simulations in steep terrain with steep groundwater gradients.

Convergence problems were encountered with MODFLOW-SURFACT and the model had to be run assuming confined conditions (while applying recharge). This illustrates the difficulty of modelling unconfined groundwater flow in steep terrain.

The model was first run with steady state conditions (using no-flow boundaries coinciding with the watershed boundaries) to evaluate regional flow and then run with transient conditions to assess groundwater contribution to streamflow under low recharge conditions.

The 3D basin-wide model was used to predict inflows to the open pit. Note that the assumption of confining conditions (required for the regional model) results in an overestimate of pit inflow estimates.

Figure 5-7: Domain of 3D basin-wide model Using MODFLOW-SURFACT), Case Study 1

5.6.1.3 Analytical Model for Pit Inflow

As fault zones were intersected during drilling in the open pit area, a Darcy-based analytical approach was used to conservatively estimate the discharge from an intersecting fault zone. For this method, flows into the pit were comprised of two components - lateral flow towards the fault and flow along the fault. These two components of flow were assumed to originate from recharge due to precipitation. Radial inflow to the pit through the surrounding lithology, which would reduce flows toward the fault, was not considered.

The estimated pit inflow along the potential fault was summed with the numerical inflow calculation and the estimate of TSF seepage contribution to determine a conservative pit inflow estimate.

5.6.1.4 TSF 2D Seepage Model

A series of two dimensional cross-sectional seepage models were created to meet the following objectives:

• Estimate seepage using a range of hydraulic conductivity values for the glacial foundation materials (sensitivity analysis)
• Identify key flow pathways
• Estimate groundwater recharge under unsaturated conditions.

The finite-element code VADOSE/W was selected for its applicability to variably saturated flow problems.

The 2D model (see Figure 5-8) incorporates hydrostratigraphic units consistent with the basin-wide model and materials representing the gradation of tailings and TSF structural elements (i.e. the core, filter, and shell of the TSF dam).

Figure 5-8: 2D cross-sectional model to assess TSF seepage (using Vadose/W), Case Study 1.

The seepage estimates determined with the 2D cross-sectional models were summed and used as recharge inputs to a 3D TSF model. This approach is conservative in a sense that the 2D model cannot account for radial flow under the TSF. Generally, flow from a relatively small pond area will radiate outwards to a larger discharge area. Summing up a series of 2D sections will cause an overlap in the recharge areas, thereby overestimating recharge.

This approach illustrates the appropriate use of conservative assumptions to simplify a 3D flow problem to a 2D flow problem.

It should be noted that the use of a conservative geometry assumption (here 2D) does not necessarily result in conservative seepage estimates. In many cases, the uncertainty in hydraulic properties of the mine waste and/or underlying aquifer material is much greater than the geometry effects. This is particularly true for variably saturated flow problems where unsaturated hydraulic properties (such as soil water characteristic curves) are required.

5.6.1.5 TSF 3D Groundwater Model

A modified version of the 3D basin-wide groundwater flow (in MODFLOW-SURFACT) was used to assess the potential for seepage from the TSF to downgradient surface water receptors.

For this purpose, the model domain was reduced in extent to exclude the steep topography west of the mine (Figure 5-9). This was done to facilitate model convergence. In addition, the tailings deposit was explicitly included in a new model layer. Furthermore, surface water features in the footprint of the TSF were removed and recharge estimates obtained from the 2D cross-sectional model were applied to the tailings beaches. The tailings pond was simulated with the MODFLOW-96 river package which permits flow between surface water features and groundwater based on a conductance assigned to the bottom of the surface water body.

Sensitivity analyses were completed using the 3D TSF model to bracket the potential range of seepage from the TSF and to estimate the relative proportion of TSF seepage reaching different surface water receptors.

Figure 5-9: Model domain for 3D TSF Model (using MODFLOW-SURFACT), Case Study 1.

5.6.2 Case Study 2: Underground Mine

Case Study 2 illustrates the use of a 3D saturated flow model to study the potential environmental impacts of an underground development. A more detailed description of this case study is provided in Appendix C.

A basin-wide fully three-dimensional saturated flow model was used to assess

• Transient dewatering during adit construction
• Transient dewatering during mine operation
• Transient simulation of post-closure rebound, flow, and transport.

5.6.2.1 Groundwater Flow Model

The finite-difference code MODFLOW was used to simulate groundwater flow. Visual MODFLOW was used as the GUI.

A 3D numerical model was essential for simulating the complex groundwater flow and discharge associated with the underground mine workings and spatial distribution of hydrogeologic units.

Although continuum models, like MODFLOW, do not address discrete bedrock fractures, the modeller noted that the applicability of an equivalent porous media approach to simulations of flow and transport in heterogeneous, fractured bedrock settings has been well proven and documented. Beyond this, the rationale for adoption of an EPM approach for this specific site was not clearly expressed in the model documentation.

Initially, a detailed 3D model had been developed which included 40 layers bounded by surface watershed divides. However, this model was found to be too unwieldy and the vertical discretization was subsequently reduced from 40 layers to 17 layers without compromising the delineation of hydrostratigraphic zones. To represent the underground mine workings, simple drains replaced the more complicated inactive zone-drain combination used in the updated model.

The updated model was also expanded to the north and south to be able to simulate potential cross-boundary flow from the adjacent watersheds to the underground mine during mine dewatering (see Figure 5-10 for updated model domain and boundary conditions). Fourteen surface water features were included as constant head and general head boundary conditions (lakes and marshes, respectively).

The revised model significantly reduced the overall complexity (and run-time) of the numerical model, yet at the same time provided more reliable modelling results. This example illustrates that a judicious use of model discretization can significantly reduce cost and time without compromising the modelling objectives.

Figure 5-10: Model domain for 3D Basin-wide Model (using MODFLOW), Case Study 2.

5.6.2.2 Solute Transport Model

The particle tracking code MODPATH and the solute transport model MT3DMS were used to simulate contaminant transport for the post-closure conditions (after mine flooding). CoCs for this project included molybdenum and arsenic which were observed in elevated concentrations in mine water discharging from an existing adit.

First, particle tracking was used to estimate the advective path and travel time of a group of particles introduced to the system (Figure 5-11). Particles were assigned to the new adit, the underground stopes, and beneath the load-out facility. Particle tracking confirmed that mine CoCs may reach groundwater discharge points several decades after closure. The particle tracking exercise helped to identify potential receptors, but was not able to provide estimates of CoC concentrations (nor travel times for reactive contaminants).

Figure 5-11: Pathline analysis using MODPATH, Case Study 2.

Following the MODPATH simulations, the solute transport module, MT3DMS, was used to simulate transport of potential CoCs in mine water to the receiving environment. No specific solute was simulated in the transport model. Instead, an arbitrary concentration of 100% was applied to all granodiorite and mine cells. The model did not simulate adsorption/absorption to rock surfaces or degradation.

This example of solute transport modelling illustrates the recommended phased approach for contaminant transport modelling. First, particle tracking analysis was completed to determine whether the CoCs present in mine water could potentially reach VECs. Second, solute transport modelling was conducted to predict the timing and magnitude of CoC concentrations at various VECs (drinking water wells).

5.6.3 Case Study 3: Groundwater Extraction Project

Case Study 3 illustrates the use of a 3D saturated flow model to study the potential environmental impacts of a large groundwater extraction project. A more detailed description of this case study is provided in Appendix C.

An existing three-dimensional groundwater flow model was used as the basis for numerical modelling. The model domain was reduced and discretization was refined to focus on the Site's surface water and groundwater systems. The three-dimensional MODFLOW code was used to simulate steady-state and transient flow.

The model domain covers approximately 80 km² of a regionally extensive aquifer. The domain boundaries were defined by large surface water features and the US-Canada border, which lies at a considerable distance to the south of the wells being assessed and therefore would not influence the flow solution (Figure 5-12). Vertically, the model was discretized into 11 layers representing the hydrostratigraphic units identified in the conceptual model. Creeks and streams within the model domain that were interpreted to be perched above the aquifer were assigned river boundary conditions, while those interpreted to intersect the aquifer were assigned drain boundaries.

The groundwater flow model was calibrated against a spatially and temporally distributed set of observations, including transient calibration to a pumping test. Predictions included an estimate of the groundwater zone of influence due to increased withdrawal rates and estimates of baseflow impacts to potentially vulnerable creeks.

This case study is a good illustration of the level of detail required to evaluate potential impacts of a large groundwater resource project in a regional aquifer which is influenced by seasonal recharge and multiple users.

Review Questions

1. What are the two most important factors in selecting the appropriate mathematical model?
   1. Modelling objectives and conceptual model.
   2. Modelling objectives and budget.
   3. Budget and conceptual model.
   4. Availability of data and modeller ability/training.
   5. Proximity of project to residential areas and valued ecosystem components (VECs)
2. The main advantages of using an analytical model to simulate groundwater flow or solute transport (as opposed to a numerical solution) are:
   1. Analytical calculations are quick and transparent.
   2. Analytical solutions are more correct and better established.
   3. Analytical solutions deal with heterogeneity and complex geometries.
   4. A and B.
   5. None of the above.
3. The main disadvantages of using a numerical model to simulate groundwater flow or solute transport are:
   1. The time and cost associated with constructing and running numerical models is significant.
   2. There is a tendency to accept numerical model simulations as the correct representation of the hydrogeologic system even though there may be considerable uncertainty.
   3. Numerical models are hard to use and are not widely accepted.
   4. A and B.
   5. A and C.
4. What flow domain approach is commonly applied to overburden material and highly weathered bedrock with high primary porosity?
   1. Discrete Fracture Network (DFN).
   2. Equivalent Porous Medium (EPM).
   3. Valued Ecosystem Component (VEC).
   4. Dual Porosity Medium (DPM).
   5. None of the above.
5. What should a model reviewer consider when determining if an appropriate mathematical model has been selected?
   1. Availability and quantity of reliable data.
   2. Dimensionality of flow and complexity of the conceptual model (geometry, heterogeneity, type of flow domain, etc.).
   3. Sensitivity of the VECs to model error from simplifications/assumptions.
   4. Time dependency and significance of transport simulation for modelling objectives.
   5. All of the above.

---

Proceed to Section 6: Numerical Model Setup