Insights from four decades of model development on the Waterloo Moraine: A review

Abstract

Modelling has played a key role in the evolution of insights on the Waterloo Moraine groundwater system as a water source for Waterloo Region. As models evolved over the last four decades, so did new insights. Starting from a simple layer-cake concept, models eventually became three-dimensional, while the focus changed from well interference to multi-layer analysis, capture zone delineation, wellhead protection areas, groundwater age and the prediction of impact due to various land uses and threats. An important insight gained was that the reliability of predictive methods depends heavily on the travel paths from source to receptor, which in turn depend on the connectivity between the various hydrogeologic units – i.e. the conceptual model. The elusive issue of predictive uncertainty in the delineation of well capture zones is approached at two scales: the local scale representing uncertainty due to heterogeneities within the individual aquifer units, and the global scale representing different conceptualizations and scenarios, where the latter is found to be the dominant uncertainty. In the case of a multi-scenario analysis involving scenarios of equal plausibility, the precautionary approach can be used to arrive at a conservative prediction. Practical applications provide insights on the dynamics of groundwater systems, the impact of changes and threats being imposed and the effectiveness of mitigation measures, as well as the time scale at which threats and mitigation measures act. These insights should help hydrogeologists to use models more effectively, and to better understand and manage predictive uncertainties that inevitably arise in dealing with complex groundwater systems such as the Waterloo Moraine.

La modélisation a joué un rôle de premier plan afin comprendre le système d’écoulement de l’eau souterrain de la moraine de Waterloo qui est utilisée comme source d’approvisionnement de la région de Waterloo. Les modèles se sont perfectionnés depuis 40 ans et ont ouvert de nouvelles perspectives. Partant du concept simple d’empilement de couche, les modèles ont progressé pour devenir tridimensionnels, et l’attention qui avait d’abord été dirigée sur l’interférence entre puits s’est portée sur l’analyse multicouche, la délimitation des zones de captage, les périmètres de protection des têtes de puits, l’âge des eaux souterraines et la prévision des effets d’occupations du sol et de menaces diverses. Un point important, il est maintenant acquis que la fiabilité des méthodes de prédiction sont fortement tributaires du trajet entre source et récepteur, lequel trajet est lui-même fonction de la connectivité entre les diverses unités hydrogéologiques – soit le modèle conceptuel. La question problématique de l’incertitude des prédictions pour délimiter les aires d’alimentation des puits est abordée à deux échelles : à l’échelle locale, qui montre l’incertitude liées aux hétérogénéités dans les unités aquifères, et à l’échelle globale qui inclut plusieurs concepts et scénarios et pour laquelle l’incertitude est plus importante. Dans le cas d’une analyse à plusieurs scénarios équiprobables, adopter une approche prudente peut permettre d’arriver à une prévision conservatrice. Les applications pratiques renseignent sur la dynamique d’écoulement des systèmes aquifères, sur l’effet des changements et sur les menaces qui s’exercent, sur l’efficacité des mesures d’atténuation, ainsi que sur l’échelle temporelle à laquelle jouent les menaces et les mesures d’atténuation. Les renseignements obtenus devraient aider les hydrogéologues à employer les modèles avec plus d’efficacité et à mieux comprendre et maîtriser l’incertitude qui entoure inévitablement les prévisions d’écoulement des systèmes aquifères complexes comme celui de la moraine de Waterloo.

Introduction

The Regional Municipality of Waterloo, with a population of about 550,000, operates 50 well fields with 126 wells supplying most of the Region’s drinking water. The majority of these wells are situated within the Waterloo Moraine. The Moraine also supports a system of ecologically sensitive streams and wetlands within the watershed of the Grand River (Veale et al. 2014, this issue). Historically, since the late 1800s, the Moraine has been a reliable and abundant source of excellent water to the residents in this thriving part of Ontario. In the early days, finding water simply meant drilling another well, without much need for systematic exploration of the source. This changed in the 1960s, when rapid industrial growth provided the motivation for a systematic exploration of the water source.

With rapid growth also came new challenges for the hydrogeologist. First, it now became of interest to know not only how much water a well would yield, but also how multiple wells would interact with each other, where the water came from, how old it is, how susceptible to contamination it could be, how urban or rural environments might affect the water source and how to find the balance between extraction and the maintenance of healthy aquatic ecosystems. Long-term sustainability became an issue of paramount importance. Thus, the need emerged for a sound understanding of the groundwater source. The traditional focus on exploitation broadened; protection of the groundwater and its quality and quantity now became important as well.

Second, there was the challenge of the extraordinary hydrostratigraphic complexities of the Waterloo Moraine, which is characterized by a mix of often-discontinuous till, sand and gravel units that would not easily fit a standard aquifer-aquitard concept. It soon became apparent that the hydrogeological tools of the 1960s were no longer up to the new challenges. New tools were needed.

This paper reviews the evolution of modelling tools that took place as hydrogeologists were confronting these challenges, and it probes into the insights that were gained in the course of this evolution.

Early models

In the 1970s, with the advance of digital computing, computer modelling emerged as a new tool for analyzing complex hydrogeological problems. The first two digital models applied to groundwater problems in Canada were the finite difference model of the Musquodoboit Harbour aquifer, Nova Scotia, by Pinder and Bredehoeft (1968), and the finite volume model of the Welland aquifer by Frind (1970), which predicted the impact of aquifer dewatering on municipal wells. Shortly thereafter, these authors joined forces to develop the first finite element groundwater model (Pinder and Frind 1972) for the analysis of complex systems. Modelling was on its way to becoming a key tool for understanding groundwater systems.

The first practical model of the Waterloo Moraine was developed in 1973 by Frind as part of a consulting project under contract with the Regional Municipality of Waterloo (International Water Supply Ltd. 1973). This was a simple two-dimensional (2D) finite element layer-cake model consisting of a confined aquifer with an overlying leaky aquitard. The objective was to predict the aquifer response at the regional scale to increases in the pumping rates of the main municipal wells, while taking into account the mutual interference of these wells. The model was calibrated to observed heads at certain observation wells, and the calibrated model was then used to predict the response of the aquifer under different pumping scenarios.

The 2D layer-cake concept had been dictated by the limited computing power available at the time. This concept required that the complex material of the Moraine could be functionally represented as a single continuous aquifer/aquitard system, which was, however, not quite in agreement with field data showing multiple layers of limited continuity. Thus, going to three dimensions became the next challenge. A simple prototype three-dimensional (3D) model was developed (Frind and Verge 1978), but experience with that model suggested that the computing power of the 1970s was not yet adequate for realistic 3D modelling of a system as complex as the Waterloo Moraine.

Quasi-3D model

The next-best approach was the quasi-3D model, in which several aquifers and intervening aquitards could be accommodated in a multiple layer-cake structure. This type of model was based on the assumptions that the contrast in hydraulic conductivity between aquifers and aquitards would be at least two orders of magnitude, allowing flow in the aquifers to be essentially horizontal and flow in the aquitards (leakage) essentially vertical, and that flow in the aquitards is at steady state. Under these assumptions, the aquifers are coupled to each other only through the leakage terms. The governing equation for each aquifer is:(1)

where h is the head, T is the transmissivity in the x- and y-directions, with T = Kb, K being the hydraulic conductivity and b the thickness; q_P represents the pumping, q_L the leakage through the aquitards and q_R the recharge. The aquifer equations are solved one by one, with the leakage flux calculated iteratively (Chorley and Frind 1978; Frind 1979).

The quest for three dimensions also spawned an interest in finding more efficient numerical approaches for achieving 3D functionality without having to deal with a fully 3D model. The governing equations of flow and contaminant transport became the focus of this interest. One novel approach was the alternating direction technique for solving the transport equation (Daus and Frind 1985; Burnett and Frind 1987). Insight into the effect of numerical errors was provided by Daus et al. (1985).

With the quasi-3D concept and its moderate demand on computing power, the modelling of multiple aquifer systems became practical. A subsequent refinement allowed for aquifer pinch-outs (Figure 1) and windows in the aquitards by connecting/disconnecting the affected layers (Blackport 1980). Provided leakage fluxes remained small relative to flow in the aquifers, the iterations generally converged quickly. The quasi-3D approach was successfully applied to the Greenbrook well field by Rudolph (1985) and Rudolph and Sudicky (1990) (Figure 1). Key data were the isopach (thickness) maps of the individual layers obtained from available field data. Where a layer thickness became zero, the layer was taken to pinch out. Hydraulic conductivities were taken to be constant for each layer. Boundary conditions were specified either in terms of hydraulic head or flux at the boundaries, based on field observations. This model can be considered the first to capture the complex multi-layered nature of the Waterloo Moraine. However, the instantaneous steady-state assumption remained a concern because, in reality, it would take a long time for steady state to be established in an aquitard after a change in pumping.

Figure 1. Quasi-three-dimensional multi-aquifer system (from Rudolph and Sudicky 1990, with permission).

The emergence of 3D modelling

In complex groundwater systems such as the Waterloo Moraine, flow is controlled not so much by the hydraulic conductivity of the permeable units, but by their continuity and interconnectedness, particularly in the vertical direction. This fundamental fact has been emphasized by Fogg (1986) and others. Bajc et al. (2014, this issue) also notes the importance of connectivity within the Waterloo Moraine system. Because quasi-3D models are limited in representing complex connections between layers, it became clear in the 1980s that, ultimately, the future of groundwater modelling would lie in 3D. It also became clear that using observed heads to define model boundaries can be problematic (heads change as pumping changes); a better way would be to use physically definable boundaries (such as a water course or a watershed symmetry boundary) wherever possible. The model should be flexible in representing irregular geometries and deformed layers with sloping interfaces. Most importantly, because the entire Moraine forms one interacting system, a direct (non-iterative) solution would be preferable.

In the late 1980s, the US Geological Survey came out with their 3D finite difference model MODFLOW (McDonald and Harbaugh 1988), which was one of the first universal 3D groundwater models. This type of model has a simple algorithm based on a rectangular discretization in the three spatial directions, which allows local refinements. A limitation is that any grid refinement introduced at a point of interest must be extended across the entire grid in all three directions, which can be inconvenient for complex applications. For this reason, MODFLOW was not considered a first choice for the study of the Waterloo Moraine with its complex stratigraphy. Nevertheless, this package has today become the most-used groundwater model worldwide. A recent version of MODFLOW is MODFLOW-2000 (Harbaugh et al. 2000).

Meanwhile, computers continued to become more powerful, and new numerical techniques emerged. One of these techniques was the preconditioned conjugate gradient (PCG) approach for solving large matrices, which made large models with many thousands and eventually millions of nodal points feasible. A particularly successful version was the PCG solver based on a block-line relaxation technique developed by Braess and König (1995), which had been specifically designed for multilayer groundwater systems with large permeability contrasts (up to five orders of magnitude) between adjacent layers.

This solver became the key component in the development of the 3D finite element model WATFLOW (Molson et al. 1995), a model specifically designed to handle multilayered systems with deformed interfaces, irregular boundaries and natural features such as streams. The model uses triangular prismatic elements (triangles in the horizontal plane); element layers are continuous, but the material properties can be varied element by element within the layer, so that different materials can be accommodated within a layer. Thus, vertical continuity in the form of aquitard windows can be represented directly through the dataset. The model is based on the governing equation:(2)

where K_ij represents the components of the 3D hydraulic conductivity tensor, S_s is the specific storage and Q_W is the pumping rate, with x_w representing the well location. Boundary conditions are specified in terms of head or flux, and the boundaries include the top and bottom of the model along with the lateral periphery. The top boundary allows specified head or specified recharge; to control mounding, a recharge spreading layer is draped over the top of the model to conduct excess recharge from low-K to high-K areas or directly to discharge boundaries (e.g. surface water). The 2D triangular grid for the 3D prism elements can be easily generated in the horizontal plane using the automatic grid generator GRIDBUILDER (McLaren 1999), which allows refining the grid in areas of interest such as near wells.

Enhancements to WATFLOW (Molson et al. 2002) included a linearized approximation for unsaturated flow processes, which allowed the top boundary to be placed at ground surface (Beckers and Frind 2000), and an automatic calibration routine (Beckers and Frind 2001), based on zoning the system and adjusting the hydraulic conductivities within each zone to obtain the best fit to observational data. The automatic calibration facility has been proven to be extremely useful in subsequent work.

The Waterloo Moraine model

Martin (1994) first applied WATFLOW to the highly complex North Waterloo aquifer system, a part of the Waterloo Moraine. The success of this model led to the creation of the fully 3D Waterloo Moraine model (Martin and Frind 1998), which subsequently became the basis for numerous studies of the Moraine groundwater system. The conceptual model was based on the complex hydrostratigraphy of the Waterloo Moraine system, characterized by the three relatively continuous till units, the Port Stanley/Tavistock, the Maryhill, and the Catfish Creek Tills, that have been identified throughout the Waterloo Moraine. (For the most recent interpretations of the Moraine geology and hydrogeology, see Bajc et al. 2014, this issue; and Blackport et al. 2014, this issue.) Between the till units, several discontinuous glaciofluvial sand and gravel units act as aquifers. The glacial sediments are underlain by bedrock. The urban areas of the Cities of Kitchener and Waterloo are for the most part situated on the eastern flank of the Moraine, which is overlain by till that helps to protect the aquifers from contamination. The core of the Moraine with its sand hills has mostly rural land use, and this is where most of the recharge to the Moraine aquifers occurs. Windows in the till along the flanks of the Moraine also allow recharge. Geochemical evidence in terms of groundwater age also supports this recharge mechanism (Stotler et al. 2014, this issue).

Unlike earlier models that focused on individual well fields and local surrounding areas, the Waterloo Moraine model of Martin and Frind (1998) was bounded by fixed boundaries that would not be affected by changes in the pumping regime, namely the Grand, the Nith and the Conestogo Rivers, and also Boomer Creek at the north and Roseville swamp at the south (Figure 2). The database available at the time (1990s) included about 4,500 borehole logs, of which 2,044, aligned along 317 local-scale cross-sections, were selected as reasonably reliable. A typical stratigraphic cross-section is shown in Figure 3. The various materials were grouped into four aquifers and four aquitards, resulting in an eight-unit conceptual model (Figure 4), which for modelling was further subdivided into 12 element layers. The lithologies identified in each borehole log were individually interpreted in terms of hydraulic conductivity K_h (horizontal) and K_v (vertical), and the effective hydraulic conductivity within each element layer was calculated as the arithmetic mean of K_h and the harmonic mean of K_v, respectively, of the material within the layer, as required to satisfy the physical laws for flow in layered systems (Bear 1972). The interface elevations and the K-values of the element layers were interpolated areally using universal kriging.

Figure 2. Waterloo Moraine model area (from Martin and Frind 1998, with permission).

Figure 3. Waterloo Moraine, typical stratigraphic cross-section (from Martin and Frind 1998, with permission). a.s.l., above sea level.

Figure 4. Waterloo Moraine conceptual model (from Martin and Frind 1998, with permission).

An important point to note from Figures 3 and 4 is that the individual layers of the conceptual model, aquifers as well as aquitards, are generally discontinuous. The discontinuities in the aquitards were interpreted as “windows” that act as major controls on the flow system. In recent work focused on the depositional history of the Waterloo Moraine, Bajc et al. (2014, this issue) have identified numerous erosional features in the till aquitards, as well as buried bedrock valleys, in the Moraine and the underlying units that account for these discontinuities. Similar evidence has also been found elsewhere; for example, van der Kamp et al. (1994) discuss a long and exceptionally narrow chloride plume in a confined aquifer in the Canadian Prairies, and Cummings et al. (2012) provide an exhaustive review of buried valley aquifers, also in the Canadian Prairies. Buried stream channels, however, are unlikely to be detected on the basis of point data consisting of single boreholes, so the original Waterloo Moraine model had to be based on the interpretation in terms of windows.

The 3D finite element grid developed for the Waterloo Moraine model is shown in Figure 5. The individual elements range in size areally from ~10 m near the wells to ~500 m in more remote areas. The grid was first generated in the horizontal plane and then projected vertically to all element interfaces, which follow the hydrostratigraphy. The model was initially calibrated by trial and error to observed hydraulic heads in the aquifer. An area-weighted recharge of 200 mm/year was initially used (later found to be too low). Water balance calculations indicated that extraction has little effect on the overall water balance. The simulated hydraulic head distribution in Aquifer 1 is shown in Figure 6.

Figure 5. Waterloo Moraine model, three-dimensional finite element grid (from Martin and Frind 1998, with permission).

Figure 6. Waterloo Moraine model, simulated hydraulic heads in Aquifer 1 (from Martin and Frind 1998, with permission). mAMSL, meters above mean sea level.

In a model study to investigate the effect of different pumping rates, Callow (1996) applied both WATFLOW and MODFLOW to the area of the Greenbrook well field (see Figure 2). This work provided insight into the importance of the boundary conditions and the advantages of the finite element method in grid refinement and in handling complex geometries. Radcliffe (2000) used the Waterloo Moraine model to study the impact of urbanization on the aquifer system within the western part of the urbanized area of the City of Waterloo, an area of rolling hills where the main aquifer is generally overlain by a protective till layer. This area is located in close proximity to the main recharge area for the regional aquifer system (Blackport et al. 2014, this issue). The model provided insight into the impact of urbanization on the flow system and the recharge reaching the aquifer.

Wellhead protection

Meanwhile, in the 1980s, the issue of groundwater protection had emerged as communities began to realize that the cost of aquifer protection is generally much less than the cost of cleaning up a contaminated aquifer. The US Environmental Protection Agency developed a formal framework for this issue by defining the concept of the Wellhead Protection Area (WHPA) (US Environmental Protection Agency 1987, 1997).

In Canada, the Region of Waterloo in the 1990s became one of the first jurisdictions to recognize the importance of groundwater protection. For the country as a whole, the defining moment for groundwater protection came in 2000, when a municipal well in Walkerton, Ontario, became contaminated by bacteria; seven people died and thousands fell ill as a result. This tragic event led to the Walkerton Inquiry (O’Connor 2002), and the Government of Ontario soon followed up on the Inquiry recommendations with detailed guidelines and regulations for source water protection for Ontario (Province of Ontario 2004, 2006).

WHPAs are normally defined on the basis of the well capture zone, which is the area within which the recharge is captured by the well. It is within this area that a potential source of contamination can become a threat to a well. A WHPA can be defined in terms of travel time from a contaminant source to the well, and it can include all or a portion of a well capture zone, depending on travel time. The delineation can be done by tracking hypothetical particles from the well in the upgradient direction through the groundwater flow field. The particle travel paths are controlled by the flow field and depend critically on the continuity between different aquifer units. One of the earliest particle-tracking routines was MODPATH (Pollock 1988), an algorithm specifically designed for rectangular block elements such as used in MODFLOW.

MODFLOW was applied to analyze groundwater flow in the Alder Creek watershed, a subwatershed of the Grand River watershed, by CH2MHill and S.S. Papadopulos & Associates (2003). In this commercial study, the surface water routing model GAWSER (Schroeter and Associates 1996) was used to generate a recharge function, which was then used as input function to MODFLOW. The model was calibrated successfully and capture zones for the municipal wells within that watershed were generated by means of particle tracking using MODPATH.

Particle tracking for irregular geometries

The first practical particle-tracking tool for finite element models such as WATFLOW was developed by Frind and Molson (2004). This code, named WATRAC, has been specifically designed for 3D heterogeneous systems with multiple layers. The algorithm is based on the semi-analytical method of Pollock (1988), and it allows particles to pass around clay lenses and change direction by refraction at material interfaces (Bear 1972). The code acts as a post-processing routine to WATFLOW. This model was applied in a commercial study by Waterloo Hydrogeologic (2000) for delineating the sensitivity zones (analogous to WHPAs) for the regional well fields within the framework of the Region’s comprehensive source water protection strategy. The present status of this strategy is discussed by Blackport and Dorfman (2014, this issue).

WATRAC is based on a steady-state flow field, with particles inserted at the well screen for each layer and tracked upgradient until they stabilize. Muhammad (2000) and Frind et al. (2002) applied an early version of WATRAC to the Greenbrook well field (see Figure 2 for location). Results for 280 years of travel time, using an improved version (Frind and Molson 2004), are shown in Figure 7. The cross-section in Figure 7 shows that the particles follow along expected paths from the well to the surface, moving preferentially through the aquifers but also passing through aquitards (note that the particle paths are off the plane of the cross-section). With some judgement, these particle tracks could be used to delineate a capture zone. We will return to capture-zone delineation in the next section after considering the uncertainty issue.

Figure 7. Greenbrook well field, particle tracks in plan view and in cross-section at 280 years (from Frind and Molson 2004).

Accounting for local-scale uncertainty: The capture probability concept

Tracking particles upgradient along a flow field does not, however, address the pervasive issue of uncertainty due to problems such as aquifer heterogeneity. Well capture zones are particularly sensitive to such uncertainties, as shown by Franke et al. (1998) and others. One way to address uncertainty due to aquifer parameters is to perform a Monte Carlo simulation, where a number of different realizations of the system are created by perturbing the uncertain parameters, and the corresponding particle tracks are then combined to create a capture zone “envelope” (van Leeuwen et al. 1998). The underlying assumption is that the uncertain parameters are normally distributed. This approach was used by Meyer et al. (2014, this issue).

Sudicky (1986) investigated uncertainty in the parameters of the porous medium in his classical field experiment at the Canadian Forces Base Borden, Ontario, which showed that seemingly homogeneous materials such as sand can be intrinsically heterogeneous. Sudicky sampled the hydraulic conductivity K of a sand unit in great detail and found log (K) to be normally distributed, allowing the material to be described in terms of its statistical parameters. On that basis, he determined that the effect of the heterogeneity on transport can be described by a macrodispersion process according to the theory of Gelhar and Axness (1983). A physical explanation of the evolution of the macrodispersion phenomenon was subsequently provided by Frind et al. (1987) through micro-scale modelling. Macrodispersion can be expressed through the standard advective-dispersive transport equation of the form:(3)

where C is the concentration, v is the pore water velocity and D_ij is the macrodispersion coefficient. Thus, the transport equation can be used to express the effect of local-scale heterogeneity and, since this heterogeneity is essentially unknown, local-scale uncertainty.

The macrodispersion approach can also be applied to capture zone delineation by applying the transport equation in a backward mode and solving for capture probability P upgradient from a well, instead of concentration C downgradient from a contaminant source. Conceptually, the advection-dispersion equation describes the physical process of advection due to a flow field, combined with dispersive spreading through the porous medium, and it can be applied to any quantity subjected to this process (e.g. solute mass, age, etc.). Capture probability has a value of 1.0 at the well screen and tends to zero far from the well in the upgradient direction. To solve for capture probability, the transport equation is switched into backward mode by reversing the sign on the advection term. The flow field is assumed to be at steady state. A formal theoretical basis is provided by the adjoint theory of Uffink (1989), Wilson and Liu (1995), Neupauer and Wilson (1999), and others.

The capture probability equation thus becomes:(4)

where τ is backward time relative to initial time of τ = 0, and v is the velocity. Because the porous medium is the same in both directions, the macrodispersion coefficient is the same as in the forward equation. The source boundary condition is P = 1 at the well. Conceptually, this approach can also be seen as a straightforward extension of upgradient particle tracking, with a dispersive component added to the advective component to account for local uncertainty. Dispersivities used in the dispersion tensor are estimated based on system scale, based on the work of Gelhar and Axness (1983), and others. The result is a plume of capture probability extending in the upgradient direction from the well. The backward transport solution is incorporated in the code WTC (Waterloo Transport Code; Molson and Frind 2005).

The capture probability plume can be interpreted as the risk level of a parcel of water being captured by the well in the presence of uncertainty. Near the well, the risk of capture is near 100%, while more distant from the well, it would gradually decline to 0%. Such an approach is conceptually more realistic than the conventional approach based on a capture zone obtained from an envelope drawn around a bundle of particle tracks, with the interpretation that the well captures all the water from inside the envelope, but none from the outside. In terms of land use, the conventional approach would call for certain restrictions on areas inside the WHPA, while no restrictions would apply to areas outside. This approach can lead to conflicts if the location of the WHPA outline is subject to uncertainty. The capture probability approach eliminates this potential problem.

The capture probability approach was applied by Frind et al. (2002) to the Greenbrook well field, using the Martin and Frind (1998) Waterloo Moraine model. Figure 8 shows the resulting probability plume at 40 and 280 years in terms of the maximum probability values over depth, as well as the cross-section at 40 years, where 280 years corresponds essentially to steady state. The cross-section in Figure 8 shows how the plume passes through a window from Aquifer 2 (the pumped aquifer) to the overlying Aquifer 1. Frind et al. (2002) extracted a 2D capture zone from the 3D probability plume by selecting the contour that would enclose an area within which the recharge would balance the pumping; this was found to be the 0.25 contour. This contour, emphasized in Figure 8, is found to be considerably wider and longer than the particle tracks (Figure 7), extending partly around a neighbouring well to the west, thus demonstrating the effect of uncertainty. The tips of the plume can be associated with a low capture probability. The capture zone represented by the 0.25 capture probability contour can be seen as being on the safe side as far as protecting the Greenbrook well field is concerned. The authors recommended that particle tracking and backward transport be used together to combine the advantages of both methods.

Figure 8. Greenbrook well field, capture probability plumes at 40 and 280 years, maximum value over depth; bottom: cross-section at 40 years (from Frind et al. 2002, with permission).

Groundwater age and life expectancy

The Clean Water Act of the Province of Ontario (Province of Ontario 2004, 2006) requires WHPAs to be delineated in terms of time of travel (TOT) of the groundwater of 2, 5 and 25 years, as well as the ultimate capture zone, allowing the characteristics of different types of contaminants to be accounted for. TOT zones are normally delineated by means of tracking particles backward for the required time. One way to account for uncertainty would be by performing a transient transport solution for capture probability and recording the plumes corresponding to the required travel times.

A further alternative approach that attaches uncertainty directly to travel time is based on the concept of groundwater age and life expectancy (Goode 1996; Bethke and Johnson 2002), which has been extended to aquifer and wellhead protection by Molson and Frind (2004, 2012). In age modelling, the mass of water being transported through the pore space by advection and dispersion is considered, with age A being a physical attribute that grows at a rate of 1 day per day, and where age is measured relative to the time of recharge at ground surface. Dispersion again expresses local-scale uncertainty, and the dispersion parameters are the same as in the standard transport equation or the probability equation. In a closed system, age will increase indefinitely, but in a groundwater system, young water continually enters via recharge and mixes with the older water, so the age distribution will eventually reach a steady state.

The equation defining this steady-state age distribution is of the same form as the standard transport equation (Equation 3), albeit with a growth term added:(5)

which is solved with a fixed source condition of A = 0 along the recharge boundary. In the same way as a particle can be tracked either forward (downgradient) or backward (upgradient), the above equation can be applied in the forward or backward direction with respect to flow. In the backward mode, instead of age A with respect to a source, the life expectancy E is solved with respect to a sink such as a well, expressing the time remaining for a particle before being captured by the sink. Life expectancy also increases with distance (upgradient) from the well. The equation is of the same form as the age equation, except for the sign on the advection term: (6)

which is now solved with a sink boundary condition of E = 0 at the well.

The steady-state life expectancy plot with respect to a given well defines all possible TOT zones for that well, again accounting for local-scale uncertainty. For any given point within a flow system, age and life expectancy will add up to the total travel time from the point of recharge to the well. Molson and Frind (2012) found the life expectancy contour for a certain value of E (from Equation 6) to be in close agreement with the P = 0.5 contour of the probability plot (from Equation 4) for the same travel time. This is demonstrated by means of the Greenbrook well field example (Figure 9), which compares the life expectancy plot at steady state (Figure 9a) to the capture probability plot at 60 years (Figure 9b), showing the 60-year life expectancy contour to agree closely with the 0.5 contour of the probability plume. The macrodispersion coefficient is the same in both solutions. Both simulations are in 3D where, for plotting in the 2D plan view, the minimum value in the vertical is chosen for life expectancy, and the maximum value over the vertical is chosen for capture probability. (In a 2D simulation, the agreement between E and P is found to be exact.) Thus, the life expectancy approach can be used to delineate TOT zones under uncertainty. The method can be seen as a straightforward extension of backward particle tracking. Groundwater age in the Waterloo Moraine from the geochemical viewpoint is discussed by Stotler et al. (2014, this issue).

Figure 9. Greenbrook well field: (a) steady-state life expectancy vs. (b) 60-year capture probability (from Molson and Frind 2012).

Determining impact: The well vulnerability concept

A capture zone or WHPA expresses whether a contaminant will reach a well within a certain time, but it cannot express the actual concentration in the well water if the contaminant does reach the well. Neither does it describe the actual impact of a potential threat to a well. Even if a contaminant were to reach a well, dilution in the aquifer and within the well could lower the concentration to below-critical levels. On the other hand, critical levels of contamination (drinking water limits) could be breached before the time indicated by TOT calculations. If concentrations were to reach critical levels, a capture zone would not say much about how long it would take to reach those levels, and if the cause were a spill, it would not be able to predict the duration of the threat.

These issues are addressed by formulating the wellhead protection problem in terms of well vulnerability (Frind et al. 2006), which expresses the actual impact on a well due to a source of contamination anywhere within the capture zone. Well vulnerability is similar to aquifer vulnerability in that both of these concepts express the potential impact at a receptor due to a threat at some source, differing only in terms of the receptor. Frind et al. (2006) use the standard transport equation in both the forward and backward mode, accounting for the well by means of a boundary condition, and applying a scale factor to make the forward and backward breakthrough curves compatible.

A more direct way is to include the well as a source/sink term in both equations (Neupauer and Wilson 2005; Rahman et al. 2006; Rahman, 2008). The forward equation is: (7)

where q₀ is the outflow rate at the well, q_I is the inflow rate at the contaminant source, C_I is the source strength and θ is the porosity. The backward equation becomes:(8)

where P is the travel time probability density function and P_I is the strength of the probability pulse at the well. The sink term in the forward equation becomes the source term in the backward equation. If the forward breakthrough curve is normalized by the injected mass and the backward travel time curve is scaled by the well pumping rate, both curves will be equivalent and of dimension [L⁻³].

A typical breakthrough curve at the well is shown in Figure 10. From the breakthrough curve, we can find (1) the time of first arrival of some critical concentration such as a drinking water limit, (2) the maximum concentration to be expected, (3) the time taken to reach the maximum (this would correspond to the travel time from a capture zone analysis) and (4) the time of exposure to above-critical concentrations.

Figure 10. Typical breakthrough curve for well vulnerability mapping (Frind et al. 2006, with permission).

The backward model applies if the potential impact is sought due to possible contaminant sources upgradient of the well, which would be of interest in the case of non-point sources (for a discussion of sources see Sousa et al. 2014, this issue). In that case, a pulse is tracked in the upgradient direction from the well and the breakthrough curve is recorded for all points of interest, usually at potential source locations at ground surface. If these points include all grid points, the results can be contoured to produce a map of well vulnerability for the given well.

Rahman (2008) and Rahman et al. (2010) applied the well vulnerability approach to one of the municipal wells within the Mannheim well field operated by the Regional Municipality of Waterloo. The original database of Martin and Frind (1998) was used, with an average net recharge of 270 mm/year and a pumping rate of 6,950 m³/day. The capture zone for this well was obtained by means of the methodology described by Frind et al. (2002). Two types of well vulnerability maps were obtained: one showing the maximum concentration to be expected in the well water with respect to the concentration at a contaminant source anywhere within the capture zone (Figure 11a), and a second one showing the time taken to reach that concentration (Figure 11b). Thus, if some point (x,y) is selected in Figure 11a, the value at that point represents the maximum concentration (relative to source concentration) expected in the well water due to a unit pulse applied at the same point, while the value in Figure 11b at the same point (x,y) represents the time taken to reach that value. Indentations in the capture zone are due to other nearby wells.

Figure 11. Well vulnerability maps for Mannheim well: (a) maximum relative concentration; (b) time taken (from Rahman et al. 2010).

These plots can be used, for example, to select optimum areas for the application of Beneficial Management Practices (BMPs). These are used in agricultural areas where nutrients such as nitrate in the well water have become or may become a problem, and they usually consist of a targeted reduction of fertilizer application (Sousa et al. 2014, this issue). The well vulnerability approach involves two steps: (1) a backward vulnerability analysis (as shown above) to map the potential response anywhere within the capture zone and (2) a forward vulnerability analysis to determine the actual response due to BMP measures within selected areas, as well as the evolution of the response in time. The complete procedure is outlined by Rahman et al. (2010).

Accounting for large-scale uncertainty: Scenario analysis

Uncertainty exists not only at the scale of local aquifer heterogeneities, but also at a larger scale of the stratigraphy and the conceptual model itself, including the boundary conditions. One notable example of such uncertainty is the presence/absence of aquitard windows, which could be in the form of erosional channels. These types of uncertainty are not amenable to stochastic analysis. However, a systematic scenario analysis (Poeter 2007) with a small number of possible scenarios can provide insight into their impact.

Sousa et al. (2010, 2013) examined the effect of large-scale uncertainty due to the spatial distribution of recharge in a scenario analysis concerned with the capture zone of the same municipal well within the Mannheim well field that was analyzed by Rahman et al. (2010) in the context of well vulnerability (see preceding section). The database, however, differed from that used by Rahman as new hydrogeologic data were added to the original database for the area in the immediate vicinity of the well.

For the large-scale uncertainty analysis, Sousa et al. (2013) developed three different but equally feasible scenarios by applying three different models to generate the spatial distribution of recharge. The three models were an integrated groundwater-surface water model (HydroGeoSphere; Therrien et al. 2005), a surface water routing model (GAWSER; Schroeter and Associates 1996) and WATFLOW. The resulting recharge distributions were fed into WATFLOW to produce the flow field corresponding to each scenario. Each of these flow fields was then calibrated against the same data, resulting in three different hydraulic conductivity distributions. Because each scenario calibrated equally well, the calibrations were non-unique, and other acceptable scenarios of this system may conceivably be found. This also means that it may not be possible to identify a “best” model on the basis of model calibration.

Each of the three scenarios produced a different capture probability plume. The plume corresponding to the WATFLOW recharge scenario (average recharge 255 mm/yr, pumping rate 6,756 m³/day) is shown in Figure 12a. The 0.5 contour of this plume can be used as a capture zone outline. This outline can be compared with the capture zone of Rahman et al. (2010) (Figure 11), which is based on a similar average recharge and the same pumping rate. The two capture zones are also similar, with differences attributable to the effect of database modifications and the different recharge distribution.

Figure 12. Capture probability maps for Mannheim well: (a) single-scenario analysis (WATFLOW), protection objective; (b) multi-scenario analysis, protection objective; (c) multi-scenario analysis, mitigation objective (from Sousa et al. 2013). The 0.5 contour has been emphasized.

Sousa et al. (2013) applied the precautionary approach to the three scenarios to address the dual objective of BMPs of groundwater protection and mitigation of contaminated areas. In order to do this, the three plumes (not shown here) were merged to produce a protection map (Figure 12b) and a mitigation map (Figure 12c), where the former used the maximum capture probabilities, and the latter the minimum capture probabilities, over all scenarios. The protection map includes all areas that might contribute to a “fail-safe” zone of uncertainty, while the mitigation map includes only areas with a high probability to contribute to a “maximum benefit” zone of confidence (Evers and Lerner 1998). In this way, the precautionary approach always leads to a conservative capture zone.

The choice of the 0.5 contour of the capture probability plume is based on the life expectancy work by Molson and Frind (2012) (see also Figure 9). Other choices are possible; for example, the 0.25 contour, as recommended by Frind et al. (2002) on the basis of the mass balance (see also Figure 8). The different contour levels on the capture probability plume represent local-scale uncertainty. In the capture probability maps of Figure 12, choosing the 0.25 contour would extend the tip of the effective capture zone, but the effect would be much less than that due to the multi-scenario analysis (see below).

The effect of the multi-scenario analysis is shown by comparing Figures 12a and 12b. Focusing again on the respective 0.5 contours, it is evident that the multi-scenario analysis resulted in a major increase of the capture zone area, an expansion much greater than that due to the locally updated database (compare Figures 11 and 12a), and also much greater than the effect of local-scale uncertainties (compare different contour levels in Figures 8 and 12). This demonstrates the overriding importance of large-scale uncertainty, which can only be investigated by means of multi-scenario analysis.

The difference between the protection objective and the mitigation objective in terms of the corresponding capture zones is shown by comparing Figures 12b and 12c. Again, taking the 0.5 contours as representative, the difference between the zone of uncertainty in the context of protection (Figure 12b) and the zone of confidence in the context of maximizing the impact of BMP measures (Figure 12c), under the selected scenarios, is huge. While this result should not be generalized, it is clear that conceptual uncertainties at the large scale can have a large or even dominating impact.

A further comparison can be made in terms of the approach used to represent local-scale uncertainty. Meyer et al. (2014, this issue) applied a Monte Carlo approach to delineate the composite capture zone for the Mannheim West well field, which contains the well analyzed above. The Monte Carlo approach is conceptually equivalent to the capture probability approach discussed above. Although the two models are not quite the same due to differences in the database and the pumping rates, both are based on a single conceptualization. Local-scale (parametric) uncertainty is accounted for in the Monte Carlo approach by random perturbations of the model parameters, and in the capture probability approach by introducing macrodispersion. Comparing Figure 12a to Figure 9 in Meyer et al. (2014) reveals that the capture zones obtained by the two approaches are remarkably similar in orientation and extent. The differences are much less than those due to the multi-scenario analysis discussed above (Figures 12a and b). This suggests that, regardless of the methodology used, local-scale uncertainties are of secondary importance. The effect due to the different methodologies should be investigated further by applying both approaches to the same database under the same scenario.

The above local- and global-scale uncertainties exist in addition to the usual uncertainty due to incomplete data, which can be remedied by collecting additional data. The Clean Water Act of Ontario (Province of Ontario 2006) addresses this uncertainty by calling for the concept of the “living model”, which is continually updated (and recalibrated) as new and better data become available. Thus, under the living model concept, the “best” model will be the most recent one. On the other hand, the multi-model analysis concept discussed above is based on different conceptualizations (for example, recharge distributions) where all are equally plausible, so a “best” model cannot be identified, but the precautionary approach can be used to arrive at a conservative result. Neither the global-scale uncertainty nor the local-scale uncertainty, as discussed in the foregoing, can be addressed by collecting additional data, at least not by practical means.

Addressing a threat: Aquifer contamination due to road salt

One of the threats to the quality of groundwater discussed by Sousa et al. (2014, this issue) is road salt. The Greenbrook well field (see Figure 2 for location) has seen chloride concentrations in the well water increase gradually in the 1990s. This well field was originally situated at the edge of the City of Kitchener, but due to increasing urbanization it is now surrounded by the City, with numerous arterial and suburban roads as well as a nearby expressway. Thus, road salt was considered to be the most likely source of the chloride.

To address this threat, the Region initiated a Road Salt Management and Chloride Reduction program to identify options for stabilizing or reducing the salt input. A key component was a chloride modelling study by Bester et al. (2006). That study utilized the transport model WTC (Molson and Frind 2005), which is fully compatible with the flow model WATFLOW. For transport modelling, a subgrid covering the Greenbrook area was cut from the Waterloo Moraine grid. Available data consisted of the record of chloride concentrations in the well water, plus estimates of salt usage per km for the various classes of roads, including snow dumps. From these data, a chloride input function was constructed corresponding to an average value of about 1400 mg/L and an average recharge of 250 mm/year. This input was then adjusted by calibrating the flux-weighted average concentrations calculated by the model to the observed chloride concentrations in the Greenbrook well water. The model was used to simulate the historical evolution of the chloride plumes emanating from the road network, giving a graphic picture of the dramatically increasing contamination of the groundwater, with plumes approaching the wells from various directions (Figures 13 and 14).

Figure 13. Greenbrook well field, chloride distributions in plan view at (a) ground surface, (b) detailed view at ground surface and (c) in pumped aquifer (from Bester et al. 2006, with permission).

Figure 14. Greenbrook well field, chloride distributions, cross-sectional view (from Bester et al. 2006, with permission). Arterial roads indicated. Cl, chloride.

The simulation made it clear that if no action were taken, the chloride concentrations in the well water would continue to increase, and even if all the salt input were to cease immediately, concentrations would continue to increase for several years before declining because of the salt already in the system. Aquitard windows were found to play a key role in the transport process. Various options for controlling the chloride input (see Bester et al. 2006) were investigated, ranging from a reduction of road salting within the entire area to total elimination of salting in certain critical areas, or elimination of salt from secondary roads. The results provided key input into the Region’s Road Salt Management Program, which was subsequently adopted.

Integrated surface-subsurface modelling

Groundwater and surface water are part of the hydrologic cycle. The exchange between these two domains is continuous and works both ways. Most conventional groundwater models place an interface between groundwater and surface water domains and treat surface water influx as a boundary condition, while most surface water models consider groundwater as simply a loss term. Some models couple a surface water module to the groundwater module to determine the exchange flux crossing the interface; this approach generally works well when groundwater flow is the dominant process. However, in situations where both groundwater and surface water are important and the two-way exchange is significant, integrating groundwater and surface water within one model is the better way.

One of the most advanced fully integrated groundwater-surface water models is HydroGeoSphere (HGS) (Therrien et al. 2005). This model uses the diffusion-wave equation to represent surface water flow, together with Manning’s equation to compute flow velocities, Richards’ equation with the van Genuchten formulation for unsaturated flow and the standard groundwater equation for flow below the water table. HGS has been applied to challenging problems ranging from flow in the Grand River watershed to the impact of glaciations on continental groundwater flow systems (Lemieux et al. 2008). An interesting feature of this type of model is that surface water bodies such as streams or wetlands need not be specified directly – the model forms these bodies automatically on the basis of the digital elevation map (DEM) provided.

The first application of an integrated model in the Waterloo Moraine area was by Jones et al. (2008), who used InHM (VanderKwaak and Logue 2001), a precursor of HydroGeoSphere, to model the Laurel Creek subwatershed in the northern part of the Waterloo Moraine. The authors give a comparison between the calculated and the actual drainage network within the Laurel Creek watershed, showing excellent agreement. They also show a near-perfect match between the computed and observed discharge hydrographs for one of the stream gauging stations in the watershed, demonstrating that baseflow as well as storm runoff are modelled with very good accuracy.

Jones et al. (2009) used HGS in a comparative study of unsaturated flow processes together with FEFLOW (Diersch 2006), a popular commercial finite element model that provides a convenient user interface. The authors applied both models to the Alder Creek subwatershed within the Waterloo Moraine. In order to achieve an accurate representation of the unsaturated zone, the model was finely discretized in the vertical, with 87 layers between ground surface and bedrock. To represent unsaturated zone processes, HGS uses the standard van Genuchten formulation, while FEFLOW uses the simpler modified van Genuchten form. A comparison of calculated heads shows that differences are small below the watertable, but much larger above the watertable. The effect of these differences on the delineation of well capture zones remains to be investigated.

Conclusions: What we learned and where we are today

By providing a rich source of data, the Waterloo Moraine has served as a living laboratory in the evolution of increasingly sophisticated modelling tools for the analysis of groundwater problems. These problems have grown in scale from the well scale to the scale of the entire Moraine system, and in complexity from the basic water supply problem to the delineation of wellhead protection areas and the assessment of well vulnerability. To meet these challenges, it became clear that the hydrogeological structure of the entire Moraine system, particularly the interconnections between the numerous different stratigraphic layers, must be understood. As understanding and insight evolved, modelling tools have also evolved from simple layer-cake models to 3D models, particle-tracking models, forward/backward transport models, and integrated groundwater-surface water models. Our review paper chronicles the evolution of these models, while the companion paper by Meyer et al. (2014, this issue) presents the current state-of-the-science for modelling the Waterloo Moraine.

A profound insight that came out of our work with the Waterloo Moraine is that uncertainty plays a key role in groundwater modelling, and that it acts at more than one scale. At the local scale, uncertainty is due to the unknown distribution of materials within an aquifer unit, while at the larger global scale, it is due to different model conceptualizations, for example, the presence/absence of windows, or the spatial distribution of recharge in the subsurface. By far the greater impact is due to the large-scale uncertainties arising from the conceptual model. These uncertainties should be understood to exist in addition to the usual uncertainties due to insufficient field data, which are remedied by collecting more data (usually constrained by the given resources) and recalibrating the model. This is a process of continual improvement, with the end result commonly accepted as the “best” model. On the other hand, conceptual uncertainties involving multiple equally plausible conceptual models will not lead to a single “best” model; however, the precautionary approach can be used to arrive at a conservative prediction.

Finally, practical applications have demonstrated that modelling is a powerful (and sometimes the only practical) tool for understanding the dynamics of groundwater systems, the impact of changes and threats being imposed and the effectiveness of mitigation measures, as well as the time scale at which threats and mitigation measures act. This understanding is particularly important for complex systems such as the Waterloo Moraine.

The above insights, particularly on uncertainty, are intended not as a discouragement but a challenge to model users. Uncertainty in its many facets cannot be eliminated, but it can be understood and managed, provided that modelling is paired with common sense. It is hoped that the insights offered in this review will help the user to find the right balance.

Dedication

This paper is dedicated to the memory of Bob Farvolden and John Nunan, pioneers of groundwater science and practice in Canada, who launched the first author on his career in groundwater modelling, and who also influenced the careers of the second, third and fourth authors.

Disclaimer

It should be noted that any well capture zones, wellhead protection areas or well vulnerability zones shown in this paper are for illustrative purposes only and do not necessarily represent the official source protection areas of the Regional Municipality of Waterloo. For official designations, please refer to the relevant documentation of the Regional Municipality of Waterloo.

Acknowledgements

This paper would not have been possible without the generations of graduate students who contributed to the evolution of groundwater modelling over the past four decades. Key parts of this review came from the work of Michelle Bester and Rengina Rahman. Hydrogeologists at the Region of Waterloo generously shared key data. The reviewers and the editors of the Journal improved the paper through insightful suggestions. The authors thank all of them. The original underlying research has been supported by the Natural Sciences and Engineering Research Council of Canada, the Province of Ontario and the Regional Municipality of Waterloo. Financial support for the production of this review has been provided by the Groundwater Geoscience Program, Geological Survey of Canada, Natural Resources Canada. The authors thank the National Ground Water Association, Elsevier Publishers, the International Association of Hydrological Sciences and Canadian Science Publishing or its licensors for permission to re-use published material.