A three-dimensional groundwater flow model of the Waterloo Moraine for water resource management

Abstract

The Waterloo Moraine has been a drinking water source for the people of Waterloo Region for over a century and, as such, it has been the subject of numerous geologic and hydrogeologic studies for over five decades. Two of the companion papers in this Special Issue describe the evolution of the hydrogeological conceptualization of the Moraine sediments and the history of groundwater modelling of the Moraine groundwater flow system, respectively. This paper builds on those findings and describes the development and calibration of a three-dimensional finite-element groundwater flow model. A key aspect in the development was the implementation of a spatial geodatabase that links the conceptual hydrogeological framework with the numerical groundwater flow model. The model was based on a detailed characterization of the groundwater and surface water systems, and calibrated to available data under average (steady-state) and variable (transient) pumping and climate conditions. Following model development and calibration, the model was used to conduct a detailed water budget and risk assessment study that compared groundwater demands to available supplies within the Central Grand River Watershed, a subwatershed of the main Grand River Watershed. Several scenarios involving future municipal water demands and potential reductions in groundwater recharge due to planned land-use development were simulated, leading to the conclusion that the projected municipal water demand to 2031 can be supplied by the existing system of wells without causing a significant reduction in groundwater discharge to ecologically sensitive streams and wetlands. The model was also applied to delineate the capture zone for a well field in the Region under conditions of uncertainty, demonstrating a methodology that could be applied to other well fields. The model provides an effective and efficient tool for Regional water managers for the long-term sustainable management of the groundwater resources of the Waterloo Moraine.

La moraine de Waterloo constitue une source d’eau potable pour la population de la région de Waterloo depuis plus d’un siècle, et elle a donc fait l’objet de nombreuses études géologiques et hydrogéologiques. Deux des articles de ce numéro spécial décrivent d’une part l’évolution de la conceptualisation hydrogéologique des sédiments de la moraine, et d’autre part l’histoire de la modélisation du système d’écoulement des eaux souterraines de la moraine. Basé sur ces résultats, le présent article décrit le développement et le calage d’un modèle tridimensionnel d’écoulement des eaux souterraines par la méthode des éléments finis. Un aspect clé de ce travail de développement a été la mise en place d’une base de données géospatiales qui établit le lien entre le cadre hydrogéologique conceptuel et le modèle numérique d’écoulement. Ce modèle a été basé sur la caractérisation détaillée des systèmes d’écoulement de l’eau souterraine et de surface, et il a été calé à partir des données disponibles dans des conditions climatiques et de pompage moyennes (stables) et variables (transitoires). Une fois le développement et le calage du modèle terminés, celui-ci a été utilisé pour une étude détaillée d’évaluation du risque et du bilan hydrique qui compare la demande en eau souterraine à l’approvisionnement disponible dans le bassin versant central de la rivière Grand, un sous-bassin du bassin versant principal de la rivière Grand. Plusieurs scénarios ont été simulés portant sur la demande future d’eau par les municipalités et de réduction potentielle de la recharge des eaux souterraines en raison notamment d’une plus grande utilisation des terres, ce qui permet de conclure que la demande d’eau municipale prévue jusqu’en 2031 peut être assurée par le réseau de puits existant, sans entraîner de réduction importante de la décharge de l’eau souterraine vers les zones humides et les cours d’eau écologiquement sensibles. Le modèle a été utilisé afin de délimiter la zone de captage du champ de puits dans la région, dans des conditions d’incertitude, ce qui démontre que la méthode serait applicable à d’autres champs de puits. Le modèle constitue, pour les gestionnaires des ressources hydriques de la région, un outil efficace et fonctionnel pour la gestion durable et à long terme des ressources en eau souterraine de la moraine de Waterloo.

Introduction

The Regional Municipality of Waterloo (the Region) lies in south-central Ontario and is home to over 550,000 residents. Approximately 75% of the Region’s municipal water demand is serviced by over 100 active groundwater supply wells located in 40 well fields, with the Grand River supplying the remaining 25% of the municipal water demand (Figure 1). The municipal water supply system serves the urban areas of Cambridge, Kitchener and Waterloo as well as portions of 16 smaller rural communities. Given the Region’s reliance on groundwater, the geology and hydrogeology of the Waterloo Moraine and its surrounding areas have been the focus of over 400 technical studies and research since the 1950s (Blackport et al. 2014, this issue). These documents provide a wealth of data and interpretations pertaining to the geology, hydrogeology and hydraulics of the aquifers and aquitards that underlie the Region. A large number of these studies were consulted during the development and calibration of the Region’s current three-dimensional groundwater flow model.

Figure 1. Region of Waterloo study area.

Frind et al. (2014, this issue) outline the history and evolution of groundwater modelling studies within the Region, including the development of models designed for capture zone delineation, well field interference and contaminant migration assessments. This paper describes a state-of-the-science groundwater flow model developed for the Region of Waterloo area, as part of a Source Protection initiative led by the provincial government, in support of Ontario’s Clean Water Act (Province of Ontario 2006). This paper outlines the approach taken to link the conceptual geologic and hydrostratigraphic models to the numerical model, and to streamline the model calibration process for current and future modellers. The application of the calibrated model for a water budget and capture zone delineation assessment in the Region is also discussed.

Conceptual hydrostratigraphic model of the Waterloo Moraine

The Waterloo Moraine is an ice-contact deposit that lies within the west-central portions of the Region of Waterloo. The Moraine rises 50 to 100 m above the surrounding till plains (Chapman and Putnam 1984) and consists of thick sequences of coarse-grained sediment with discontinuous interbeds of silt, clay and fine-grained till (Farvolden et al. 1987; Karrow 1993). The stratified sediment was deposited when ice lobes of the Laurentide Ice Sheet advanced independently out of the Great Lakes basins into the Waterloo area during the Late Wisconsinan glaciation (Karrow and Paloschi 1996). The ice lobes converged, forming a complex interlobate moraine consisting of sediment deposited in subglacial, glaciolacustrine (quiet water basinal, subaquatic fan and deltaic) and glaciofluvial (braided stream and tunnel valley) environments (Karrow 1988; Karrow and Paloschi 1996; Bajc and Shirota 2007; Russell et al. 2007). These varied environments contributed to the complex stratigraphy present within the Moraine.

The first three-dimensional geologic/hydrostratigraphic conceptual model of Waterloo Region was created by Bajc and Shirota (2007), using overburden sediment exposures, stratigraphic descriptions from continuously cored boreholes and water well records, seismic profiles, mapped surficial sediments and subsurface interpretations of other researchers. Data were compiled and interpreted to create an 18-layer hydrostratigraphic model that spans the entire Region. Geologic rules were developed to prevent the geologic model layers from crossing, or extending above ground surface or into the bedrock. A recent update of the hydrostratigraphic model is presented as a contribution to this Special Journal Issue by Bajc et al. (2014). The update represents a refinement of the structure as well as the outline of the Waterloo Moraine, but it does not affect its major architectural components. The three-dimensional flow model developed in this paper is based on the original interpretation by Bajc and Shirota (2007).

Numerical model development

The Region’s three-dimensional groundwater flow model discussed in this paper was built using the FEFLOW finite element software program (Diersch 2006). The model domain includes the entire Region of Waterloo area, plus an additional area that was set at a distance such that the external model boundary conditions would not influence the model predictions (Figure 1). At its widest, the model is approximately 70 km wide from west to east, and 65 km from north to south, giving the model a total area of approximately 2876 km².

Model design

The finite element mesh in the FEFLOW model contains approximately 192,000 triangular prismatic elements in each model layer, and approximately 96,000 nodes on each model slice (a nodal layer between two element layers). The mesh was refined in areas where an enhanced definition of groundwater flow and the potentiometric surface was desirable. Element lengths of 5 m were applied around municipal and non-municipal pumping wells that extract relatively large volumes of groundwater, and elements ranging in size from 5 to 25 m were applied along surface water features of interest such as streams and creeks (Figure 2). Element lengths varying from 100 to 400 m were applied outside the central portion of the model domain.

Figure 2. Finite element model mesh refinement (numerical model Layer 1).

Model layers

In addition to the horizontal discretization, the model domain was also discretized vertically into 21 numerical FEFLOW model layers that represent variable hydrogeologic conditions with depth. The overburden units are represented in the FEFLOW groundwater flow model with 13 model layers, and they are underlain by eight Paleozoic bedrock model layers (Figure 3). As outlined previously, geoscientists at the Ontario Geological Survey subdivided the overburden beneath the Region into 18 stratigraphic layers (Bajc and Shirota 2007). Several of these overburden layers have limited thickness and spatial extent within the Waterloo Moraine area, so several of the 18 overburden hydrostratigraphic units were combined to create 12 conceptual overburden hydrostratigraphic layers (Figure 3). Those layers were then migrated into the FEFLOW groundwater flow model; however, one conceptual hydrostratigraphic layer (Middle Waterloo Moraine Stratified Sediments; Bajc and Shirota 2007) was split into two numerical model layers. The elevations at the top and bottom of the 13 overburden model layers were updated within the urban municipal well field areas through the generation and interpretation of over 100 local-scale cross-sections (Stantec Consulting Ltd. 2009, 2012a, 2012b, 2012c; Golder Associates Ltd. 2011a, 2011b, 2011c; AquaResource Inc. 2012; Blackport Hydrogeology Inc. 2012a, 2012b). The cross-sections were interpreted using hydraulic head data, hydraulic test data, lithologic descriptions in borehole logs, downhole geophysical profiles, groundwater chemistry data, understanding of the glacial history of the area and surficial geophysical datasets. The regional-scale stratigraphic layers proposed by Bajc and Shirota (2007) were refined in the groundwater model using hydrogeology-based interpretations to ensure local-scale hydraulic connections between production aquifers, and between the production aquifers and nearby surface water features. This step was essential as some units predicted to be aquitards in the regional-scale stratigraphic model were found to be municipal water supply aquifers. As the stratigraphic layers were combined to create the FEFLOW numerical model layers, a link between the overburden conceptual geologic strata and the numerical model layers was maintained.

Figure 3. Conceptual hydrostratigraphic and groundwater flow model layers (adapted from Bajc and Shirota 2007; reproduced and adapted with permission, Queen’s Printer for Ontario).

In addition to the overburden model layers, seven bedrock layers were created for the numerical groundwater flow model through the generation and interpretation of cross-sections, primarily east of the Waterloo Moraine, where overburden thins and bedrock is the dominant source of municipal water supply. Model layers were selected on cross-sections using higher quality borehole logs, downhole geophysical profiles, groundwater chemistry data, interpretations of hydraulic test data and an understanding of the Paleozoic bedrock of the area as interpreted by Brunton (2008, 2009).

The development and interpretation of the local-scale overburden and bedrock cross-sections within the urban area of the Region represented a significant step forward in the groundwater flow modelling of the Waterloo Moraine.

Model boundary conditions: Surface water features

Potential interactions between groundwater and surface water features such as rivers, streams, ponds, reservoirs and wetlands can be simulated in groundwater flow models using fixed-head boundary conditions, with variable resistance between the surface water body and the underlying aquifer. The application of a fixed-head boundary condition to represent a surface water feature is a practical approach that is frequently applied in regional-scale groundwater flow models. The inherent assumptions with this approach include:

1. fluctuations in stream or river stage are relatively short in duration and can be represented with a long-term average value,

2. seasonal fluctuations in stream or river stage do not substantially affect the hydraulic gradient to/from the underlying aquifer, and

3. stress changes simulated in the model, such as groundwater pumping, will not cause the surface water body to be fully drained.

In southern Ontario, significant changes in stream stage (i.e. > 0.5 m) often occur for short durations, including the spring melt (1 to 2 weeks), during storm events (hours or days) and extreme dry periods (usually limited to a few weeks). Outside of these periods, the fluctuations in the stream stage are generally less than 0.5 m, a value that is small as compared to the head difference within the primary aquifer (i.e. > 5 m for areas of perennial discharge). As a result, the constant stage assumption is considered reasonable and considered to have a minor impact on simulated gradients between the aquifer and the surface water features.

Flow resistance between the aquifer and the stream is accommodated in the groundwater flow model by applying spatially variable discrete hydraulic conductivity zones around the surface water bodies. Using this approach, the applied boundary condition behaves as a head-dependent flux boundary. Application of these surface water boundary conditions allows groundwater to discharge into surface water features, or surface water to recharge the underlying aquifer, depending on the model-simulated head difference between the aquifer and the surface water feature.

Care must be applied when using specified heads to represent surface water features to ensure that any volume of water supplied via the boundary condition is less than the volume that would typically flow in the stream. This is accomplished through manual assessment of the local mass balance during the calibration and simulation process and limiting flow from surface water bodies where necessary (flow-constrained head boundary conditions). The only way to avert the need for such local mass balance checks would be to apply a fully integrated groundwater-surface water model; however, this was not practical for this application as such a numerically-demanding approach would limit the ability to calibrate the model to local hydrogeologic conditions.

Using this approach, seasonal changes in groundwater/surface water interaction, as well as those due to stress (i.e. recharge or pumping), are readily simulated. For example, Figure 4 illustrates the resultant model simulated model simulated spatial and temporal variability of hydraulic gradients along a selected stream reach of Laurel Creek, one of the tributaries of the Grand River.

Figure 4. Spatial and time-varying aquifer-stream hydraulic gradient and baseflow discharge along Laurel Creek.

Efforts were undertaken to ensure that surface water features of interest, especially those located in the urban areas near municipal wells, were well represented within the groundwater flow model. Surface water features with perennial groundwater discharge (e.g. baseflow) that lie within 1 km of a municipal water supply well were simulated in the groundwater model using specified head boundary conditions. Similarly, rivers observed to host coldwater fish communities reliant on groundwater discharge for their survival were also simulated in the model in the same manner. The stream bed elevation was estimated using a high-quality digital elevation model of the surface topography that was inspected in a geographic information system (GIS) to ensure the river stage specified in the model decreased monotonically in the downstream direction for all rivers. Where required, the streambed elevation was corrected to ensure the stream declined in a downstream direction, so adjacent boundary conditions would not erroneously supply and/or remove large volumes of water from the groundwater flow system. This cycling of water between adjacent boundary conditions can lead to model instabilities and water budget errors and, therefore, it was important to use care when assigning boundary condition values.

Model boundary conditions: Recharge

Groundwater recharge refers to the amount of water that infiltrates through the unsaturated zone to the underlying groundwater table. The rate of groundwater recharge is dependent on a number of factors, including climate, land use and vegetation, surficial soil type (geology), physiography and ground surface topography.

Recharge for the Regional Model was estimated using output from the physically based Guelph All-Weather Sequential-Events Runoff (GAWSER; Schroeter and Associates 1996) continuous streamflow-generation model. The GAWSER model utilizes Quaternary geology, land cover and vegetation to subdivide the model area into hydrologic response units that predict how that land unit will respond to a precipitation event. Precipitation, estimated using input from historic and current climate records, was partitioned into three major hydrologic components: evapotranspiration, runoff and recharge. The model was calibrated by comparing the GAWSER model-simulated hydrographs to observed streamflow at various gauge locations within the Region.

The average annual recharge applied within the groundwater flow model domain (Figure 1) was 165 mm/year and ranged from a low of 0 mm/year on groundwater discharge sites such as wetlands to a high of over 500 mm/year in areas where coarse-grained gravels were mapped at surface. Recharge rates of over 700 mm/year were assigned in the model for active, or former, sand and gravel operations where vegetation, runoff and evaporation are limited.

Model properties: Hydraulic conductivity values and zonation

In standard modelling practice, hydraulic conductivity values are assigned to cells or elements of the numerical model, or according to layers or zones of uniform hydraulic conductivity. Initial values are estimated and subsequently adjusted by calibration to observational data. For the Waterloo Moraine flow model, initial hydraulic conductivity estimates were based on the conceptual hydrostratigraphic model (Figure 3).

The in-depth understanding developed in the study of the regional-scale depositional environments played a major role in the assignment of the material properties to the model. Interpreted bedding characteristics, grain size and textural characteristics (Bajc and Shirota 2007; Karrow 1993) that may impact the vertical or horizontal flow of groundwater through each respective unit were recorded on a layer-by-layer basis. The geologic and hydrostratigraphic notes were stored spatially in a geodatabase, and were used by groundwater flow modellers during the model calibration process. For example, the Lower Maryhill Till has been described as a fine-grained clay diamict unit with minor clasts, deposited in a glaciolacustrine environment (Bajc and Shirota 2007). As the matrix of the till is clay-rich with few clasts, the initial hydraulic conductivity value of the till was interpreted to be on the order of 1 × 10⁻⁹ m/sec. Similarly, a horizontal to vertical anisotropy ratio ranging from 2 to 10 was interpreted to represent minor interbeds of silt or sand that may exist within the unit. All of this information was stored spatially on a layer-by-layer basis to aid in the model calibration process. Figure 5 illustrates the initial hydraulic conductivity zones applied in Layer 6 of the numerical groundwater flow model (Middle Waterloo Moraine Stratified Sediments; see Figure 3).

Figure 5. Initial (pre-calibrated) hydraulic conductivity zones in model Layer 6. Well fields also shown.

The hydraulic conductivity zones and values noted above were subsequently dissected into smaller zones, or the values were refined using interpreted results from aquifer tests conducted on the Region’s overburden and bedrock aquifers. Aquifer test data were not re-interpreted; rather, the interpreted hydraulic conductivity, transmissivity and storage estimates for the aquifers were organized spatially in the geodatabase for future reference.

Interpreted aquifer test data were used to help parameterize hydraulic conductivity zones, as they capture hydraulic responses of the hydrogeologic system at small and larger scales, depending on the scale of the test. Pumping or shutdown test results were examined to assess potential hydraulic interconnections within, and between, production aquifers, and between the production aquifers and nearby surface water features. Hydraulic conductivity zones were subdivided in areas with interpreted aquifer test data, based on observed water level responses in monitoring wells screened in different aquifer units, or in surface water features. Transmissivity and hydraulic conductivity estimates, interpreted from pumping tests, were used to guide the assignment of hydraulic conductivity zone(s) in the aquifer(s) near the municipal production wells; this information was also stored on a layer-by-layer basis in the project geodatabase.

A bulk hydraulic conductivity analysis was undertaken that consisted of mapping borehole lithology (with associated hydraulic conductivity values related to reported geologic material type) to numerical model layers as described in Martin and Frind (1998). Horizontal and vertical hydraulic conductivity estimates were calculated based on parallel and series flow according to the materials reported in each model layer (Freeze and Cherry 1979). Using these calculated point estimates of hydraulic conductivity values, zones were delineated using trends between neighboring point estimates. This approach yielded results similar to and consistent with the estimates derived from regional-scale depositional understanding and aquifer tests; however, the bulk hydraulic conductivity analysis incorporated greater heterogeneity into the hydraulic conductivity zones.

Results of the bulk hydraulic conductivity analysis were also recorded spatially for each layer in annotated GIS files that were linked together in the spatial geodatabase. This geodatabase houses the conceptual hydraulic conductivity estimates, the interpreted results of hydraulic testing and the results of the bulk hydraulic conductivity analysis, amongst other information. The spatial geodatabase, outlined in detail below, was instrumental in calibrating the groundwater flow model and will also be essential in the efficient operation of the model by future users.

Model calibration approach and the geodatabase

The purpose of model calibration is to develop a numerical representation of the physical system that realistically reproduces observed conditions. Numerical groundwater flow models are typically calibrated by systematically and iteratively adjusting model input parameters and boundary conditions to determine the best match (within an acceptable margin of error) between the model-predicted results and field observations. The model’s ability to represent observed conditions is analyzed qualitatively to assess trends in water levels and distribution of groundwater discharge, and quantitatively to determine statistical measures of calibration. Quantitative evaluations are undertaken for both long-term average and variable pumping (transient) conditions.

As outlined previously, the initial hydraulic conductivity estimates and their spatial distributions were based on the geologic understanding of the depositional environments (Bajc and Shirota 2007; Brunton 2008, 2009). The conductivity values and zones were refined using interpreted aquifer test results, a bulk hydraulic conductivity analysis and understanding developed during the model calibration process. This dataset was also stored in the spatial geodatabase for reference by current or future hydrogeologists or engineers. Figure 6 illustrates the calibrated hydraulic conductivity zones for numerical model Layer 6 (Middle Waterloo Moraine Stratified Sediments).

Figure 6. Calibrated hydraulic conductivity zones in model Layer 6.

As part of the calibration process, the hydraulic conductivity zones and values were revised using information stored within the geodatabase. If a hydraulic conductivity value within a zone differed from the surrounding zones, or the conceptual hydrostratigraphic model, the rationale for the updated zone or value was added to the geodatabase. This was done to ensure the conceptual model is complete, and to document changes made to the model layers so reviewers or future modellers are aware of how or why zones that may appear contrary to the conceptual model were applied.

One example where a hydraulic conductivity zone in a given layer may differ from the conceptual understanding may be the presence of an erosional feature through a regional aquitard unit. Erosional features may be represented in the model using hydraulic conductivity zones with values representative of coarse-grained deposits (e.g. sand or gravel) within an aquitard layer. Within the database, the information documenting the support for such a change from the conceptual stratigraphic layer identified by Bajc and Shirota (2007) may include the borehole identification number and a geologic description that suggest the aquitard is absent, or reference to an aquifer test result that suggests hydrogeologic conditions differ from the conceptual conditions.

Tracking changes made during the calibration process within the geodatabase is a step forward in the evolution of the modelling of the Waterloo Moraine. The geodatabase allows anyone interested in the conceptual understanding of the area to see the technical basis for the hydrogeological interpretations that were implemented in the groundwater flow model. If additional data, such as a continuously cored borehole, additional geophysical data or hydraulic test data, are collected and used to refine the conceptual hydrostratigraphic model layers, the geodatabase and numerical model layers can be updated to reflect the change.

Steady-state calibration

The groundwater flow model was calibrated to long-term average conditions (steady-state) by reducing the discrepancies between the observed and model-predicted water levels within a reasonable margin of error. The steady-state calibration focused on simulating pumping conditions and average annual water levels reported for monitoring wells and other wells for the 2003 and 2006 period.

The groundwater flow system within the urban area in the Region is constantly changing as municipal wells are turned on and off to meet demand, as new wells come online and as older wells are decommissioned. On a broader scale, the total average annual groundwater demand in the urban well fields of the Kitchener-Waterloo area has been fairly consistent since the year 2005, which has allowed the system to reach a quasi-steady-state condition on an average annual basis.

Model boundary conditions and hydraulic conductivity zones and values were manually and iteratively adjusted through the calibration process before arriving at a set of values for all layers and units that produced a reasonable simulation of groundwater levels, as compared to the observed field data. In the absence of local knowledge, information from previous modelling efforts or field-based studies, literature values (Freeze and Cherry 1979; Anderson and Woessner 2002) were consulted to ensure the calibrated hydraulic conductivity values were reasonable.

As the municipal pumping rates typically vary from month to month and year to year, the associated observed water level data was divided into high, medium, medium-low and low-quality categories for model calibration. The calibration targets were categorized based on when the observed water level was collected and whether the well was an observation well, domestic well or other well. For example, water level measurements collected in observation wells with data loggers during the 2003 model calibration year were classified as high-quality targets, whereas water level measurements collected in domestic wells decades before the calibration period were considered low-quality targets.

Table 1 outlines the calibration statistics within the urban areas of the cities of Kitchener, Waterloo and Cambridge. Within each of the well field areas, the absolute mean residual was less than 4 m for the high- and medium-quality targets, and less than 5m for the remaining targets, suggesting a good match between the model-simulated water levels and observed water levels (simulating average annual municipal demands; Matrix Solutions Inc. and S.S. Papadopulos and Associates 2012). The smallest residuals lay within the urban areas, reflecting the additional hydrostratigraphic characterization, the presence of higher quality water level data and the level of effort applied during model calibration. The normalized root mean squared residuals for all quality groups in the urban portion of the groundwater flow model varied from 4.2 to 5.5%, suggesting a good match to the data. The global mass balance error was less than 0.01%, suggesting the model reached a stable and converged solution.

Table 1. Regional model calibration summary statistics.

Quality	Count	Mean residual (m)	Absolute mean residual (m)	Root mean squared residual (m)	Normalized root mean squared residual (%)
High quality	104	1.0	3.9	5.2	4.6
Medium quality	422	1.9	3.6	4.9	4.0
Medium–low quality	127	0.3	4.1	5.0	4.2

The calibration is also presented visually as a scatterplot of model-simulated and observed water levels on Figure 7. A perfect fit line (solid diagonal) and 10-m offset lines on either side of the perfect fit line (dashed diagonals) illustrate the degree of calibration to the observed values. The model replicates the observed conditions throughout the urban area. Some areas of the model domain contain perched groundwater tables that are not currently fully replicated in all instances by the model; however, additional calibration may occur in the future to better simulate these conditions.

Figure 7. Scatterplot illustrating the steady-state head calibration.

Transient calibration

Following the preliminary steady-state calibration, the model was calibrated to several transient well pumping tests or shutdown tests to increase the confidence in the model’s ability to replicate the real world response to groundwater pumping within the urban well fields. The transient calibration involved fitting the simulated drawdown and recovery at pumping wells and monitoring wells to water level responses observed during the well pumping and shutdown tests. The transient calibration aimed to match observed changes in water levels in monitoring wells that exhibited significant responses to pumping.

Transient calibration was conducted in each of the urban well field areas within the Region including the Erb Street well field. All of the municipal wells in the Erb Street well field (Figure 5) were shut down for 15 days in 2009 for well maintenance, and water level responses in nearby observation wells were monitored before and after the shutdown. Average monthly municipal pumping rates were applied in the model on the remaining municipal supply wells and other model input parameters and boundary conditions, such as recharge, were held constant during the transient calibration.

Manual adjustments to the hydraulic conductivity and specific storage values and distributions were applied until a suitable match was achieved between the simulated and observed water level changes during the well field shutdown period. Specifically, adjustments made to the hydraulic parameters in the Erb Street well field area through the transient calibration included increasing the storage values to 1.5 × 10⁻³ m⁻¹ in the production aquifer (Middle Waterloo Moraine Stratified Sediments; Figure 3). This additional storage in the unsaturated portion of the aquifer was estimated through local well field testing (Stantec Consulting Ltd. 2009). Increasing the storage values improved the fit between the observed and model-simulated water level recovery throughout the test. Other refinements made to the conductivity and storage values are outlined in a report by Matrix Solutions Inc. and S.S. Papadopulos and Associates (2012).

Figure 8 illustrates the model-simulated and observed water level recovery at three observation wells in the Erb Street well field that were monitored during the well field shutdown. The three observation wells are located less than 500 m from one of the three municipal wells and they illustrate the water level recovery following the shutdown. Water levels were not collected in the production wells during the test due to well rehabilitation; however, approximately 2 m of recovery was noted at each of the wells in mid-February, and the magnitude of this recovery was replicated at the production wells in the model.

Figure 8. Erb Street well field transient calibration example.

Refinements made to the hydraulic conductivity values and distributions during the transient calibration were transferred back to the steady-state model to ensure the steady-state calibration was maintained.

A successful calibration to field-measured conditions, while retaining consistency with the hydrogeologic conceptualization, provides a level of confidence that the model can be applied in a predictive capacity to examine potential responses to stressors, such as changes in groundwater pumping rates or groundwater recharge, or climatic variability.

Model application

Following the steady-state and transient calibrations, the groundwater flow model created for the Region can be used as a tool to guide water management decisions. The numerical model will continue to be refined and updated as borehole and water level data are collected that refine or update the existing hydrostratigraphic conceptual model. The following sections outline two projects that applied the calibrated groundwater flow model described above.

Tier Three Water Budget and Local Area Risk Assessment study

The Region of Waterloo Tier Three Water Budget and Local Area Risk Assessment (Tier Three Assessment) was initiated in 2008 to identify water quantity users and activities that may threaten the Region’s ability to supply drinking water to its residents at its future estimated pumping rates. This study was one of several technical studies conducted across Ontario to safeguard the quantity and quality of municipal drinking water sources.

The detailed water budget conducted as part of the Tier Three Assessment was completed on the Central Grand River Watershed, which contains many of the municipal wells within the urban areas of the Region (Figure 1). Three components of the water budget were quantified using the groundwater flow model: water supply, water reserve and water demand as prescribed in the technical rules and guidance for the Clean Water Act (Province of Ontario 2006; Ontario Ministry of the Environment 2009). The water supply component was based on summing the annual groundwater recharge and model-simulated cross-boundary flows, the groundwater reserve component was based on model-simulated net groundwater discharge to surface water features and water demand was the total estimated municipal and non-municipal groundwater extraction. The percent water demand was calculated as water demand divided by the difference between the water supply and water reserve.

Overall, the average annual groundwater recharge (between 1980 and 1999) for the Watershed was approximately 200 mm/year or 3560 L/s. Net cross-boundary flows were estimated to be approximately 70 L/s, leading to a total groundwater supply of approximately 3630 L/s. The groundwater reserve and total municipal and non-municipal groundwater demand within the Watershed were estimated to be 210 and 1620 L/s, respectively. Therefore, the percent water demand was estimated to be approximately 50%, or half of the water flowing through the subsurface in this area. This estimate indicates water stored within the aquifer system is not being depleted, but the municipal and non-municipal water takers within the Watershed are capturing and utilizing a large proportion of the annual supply flowing through the subsurface.

Following the water budget, the groundwater model was applied to run a set of hypothetical water quantity risk assessment scenarios. The scenarios predicted the potential impacts associated with groundwater recharge reduction due to future land use development, and future (2031) projected municipal water demands, under average annual and drought conditions (Matrix Solutions Inc. and S.S. Papadopulos and Associates 2014). While pumping tests have been conducted on an individual well or well field scale, this assessment simulated the cumulative impacts associated with future municipal demand and land-use development. Impacts were assessed based on the model-simulated change in water levels in the municipal pumping wells and the simulated change in groundwater discharge to ecologically sensitive surface water features (as per Technical Rules, Ontario Ministry of the Environment 2009). Specifically, model-simulated water levels within the municipal aquifers at each of the municipal wells were tabulated for each risk assessment scenario alongside model-simulated changes in groundwater discharge to rivers and creeks hosting coldwater fish communities. The results were used to assess the hydrogeologic and environmental sustainability of the municipal water supplies within the Region. The Tier Three Assessment aimed to identify if the future municipal pumping rates may cause pumped water levels in the municipal wells to fall below safe operating levels, or cause unacceptable groundwater discharge reductions to streams hosting coldwater fish communities (Ontario Ministry of the Environment 2009).

The risk assessment scenarios predicted the changes in water levels and baseflow under average annual and drought climatic conditions. To predict the potential impact due to drought, the model was run transiently using monthly groundwater recharge rates derived from the calibrated GAWSER hydrologic model using a 45-year (1960 to 2005) historic climate dataset.

Results of the Tier Three Assessment water budget and risk assessment scenarios suggest the forecasted municipal water demand to the year 2031 can be met with existing groundwater wells and surface water intakes. The estimated future takings are not expected to lead to pumped water levels at any municipal wells that fall below safe operating thresholds. The impacts of increased municipal demand on sensitive surface water features, such as rivers and streams hosting coldwater fisheries, is predicted to be less than 10% of the current estimated baseflow (Matrix Solutions Inc. and S.S. Papadopulos and Associates 2014); however, if planned land development proceeds without mitigating best management practices, groundwater discharge to some sensitive coldwater streams may be reduced by more than 10% of the current measured baseflow values (Matrix Solutions Inc. and S.S. Papadopulos and Associates 2014). Municipal governments within the Region require the maintenance or enhancement of groundwater recharge for all development applications, so these predicted impacts on surface water features represent a worst-case scenario, and are greater than one would expect.

Capture zone delineation

The calibrated groundwater flow model was tested for use in delineating steady-state capture zones for four wells within the Mannheim West well field (see Figure 5 for location). The combined pumping rate of the four wells is 143 L/s. In order to account for uncertainty, the software program PEST (Parameter Estimation; Doherty 2007) was used to generate 100 statistically calibrated models with unique realizations of hydraulic conductivity values. A log-normal hydraulic conductivity distribution was specified and the values were set to vary within one order of magnitude higher and lower than the calibrated hydraulic conductivity values for the original model. This constraint preserved the contrast between aquifers and aquitards and ensured that each of the 100 realizations created would be consistent with the conceptual hydrostratigraphic model, and would produce capture zones that were within the range of uncertainty associated with estimating hydraulic conductivity values. Similar stochastic approaches were used by van Leeuwen et al. (1998), Vassolo et al. (1998), Stauffer et al. (2005) and others to estimate the uncertainty in capture zone size and extent. Sousa (2013) and Sousa et al. (2013) used an equivalent approach based on a backward solution of the transport equation, where parametric uncertainty is represented by means of macrodispersion.

Each of the 100 realizations was considered statistically equivalent to the original model and therefore each predicts an equally plausible solution. Fictional particles of water were tracked backward in time from the well screen through the municipal production aquifer to the recharge area, or to another boundary such as a stream or river. This backward particle tracking was conducted in each of the 100 realizations, and the advective particle traces for each of the 100 model realizations were processed through a GIS interface to create composite capture zones that showed the range of “most common” to “least common” areas of capture for a well field under a given pumping regime.

Figure 9 shows the resulting composite capture zone for the four selected wells. The red-coloured area in the central portion of the capture zone illustrates areas where over 90 of the 100 realizations contained particle tracks that originated at one of the wells in the municipal well field and tracked backward to the recharge area. The lighter blue-grey and white portions of Figure 9 illustrate the areas where less than five (white) or 15 (blue-grey) of the 100 realizations contained particles that originated at one of the four municipal wells.

Figure 9. Composite capture zone for the Mannheim West well field.

The capture zone delineated using this methodology was an academic exercise that was not adopted for use by the Region for their capture zone or wellhead protection area delineation; however, the results provide insight into the potential effect of parametric uncertainty.

In standard modelling practice, one model or one realization is often used to delineate municipal well capture zones. The approach used in this study addresses the probability that water entering the ground within that area will travel through the subsurface to the well screen under conditions of uncertainty in the distribution of hydraulic conductivity. This probabilistic approach helps to reduce subjectivity and increase transparency in capture zone delineation.

While the probability of capture in this example was based on a reasonable number of parameter distribution realizations, the realizations were still constrained to a single conceptual model. Alternative conceptualizations or model realizations that simulate windows in aquitards overlying production aquifers, changes in the model layer structure, or recharge could also be assessed to provide additional insight into the impact of structural uncertainty on the capture zone size or shape. Such alternative conceptualizations may to lead to more significant variations in the capture zone extent than the Monte Carlo method described above with 100 realizations. As such, a combination of methods could be employed to produce well field capture zones that encompass a range of structural, conceptual and parametric uncertainty.

The effect of alternative conceptualizations is demonstrated by Frind et al. (2014, this issue), who delineated a hypothetical capture zone for one of the wells in the Mannheim West well field analyzed above. The groundwater flow model developed in that study used pumping rates, hydraulic conductivity distributions and model layers that differed from those used in the above Monte Carlo analysis, yet produced a comparable capture zone under a single-scenario analysis. However, under a multi-scenario analysis focused on the spatial recharge distribution, the capture zone changed in major ways, demonstrating the profound effect of large-scale or conceptual model uncertainties. Capture zones cannot be tested in the field and, as such, the methodologies employed in this study, or in those of other researchers (e.g. Frind et al. 2006; Sousa et al. 2013, Frind et al. 2014, this issue) can be used to identify the most probable areas of capture, which in turn can help guide water management and protection decisions.

Conclusion

A three-dimensional finite-element groundwater flow model was developed to represent the complex hydrogeological system within the Waterloo Moraine and surrounding areas. This model was used to complete a detailed water budget assessment that examined the long-term sustainability of municipal groundwater resources in the Region. The model could be used for further hydrogeological studies and to guide various water management initiatives within the area.

Several sources of geologic and hydrogeologic data were used to develop the model structure, input parameter and boundary conditions, and many of these data sources were summarized and included in a spatial geodatabase that currently houses the conceptual geologic and hydrogeologic models, and provides a link between the complex geologic model and the groundwater flow model.

The calibrated groundwater flow model provides a management tool for staff at the Region of Waterloo. It enables the municipality to continue to measure and effectively balance the potential cumulative impacts of future water demands and land-use development with preservation of the water resources and natural ecosystem functions under current and future potential climatic conditions.

Disclaimer

The well field capture zones and wellhead protection areas shown in this paper are for illustrative purposes only and do not represent the official source protection areas of the Region. For official designations, please refer to the relevant documentation of the Regional Municipality of Waterloo.

Acknowledgements

Financial support was provided by the Groundwater Geoscience Program, Geological Survey of Canada, Natural Resources Canada. The authors thank the Ontario Geological Survey for permission to re-use published material. The authors would also like to thank Gaelen Merritt, Laura Weaver and Doug Anderson of Matrix Solutions Inc. for their roles in helping refine, calibrate and apply the groundwater flow model. Our thanks are also extended to Sam Bellamy (Matrix Solutions Inc.) for his work in the GAWSER hydrologic model calibration and discussion of the recharge distribution. Many thanks are extended to Eric Hodgins and Richard Wootton at the Region of Waterloo for their technical input during the Tier Three Water Budget Assessment. Drs. Emil Frind, John Molson, Ken Howard and Alfonso Rivera are gratefully acknowledged for their comments on earlier manuscripts, which greatly clarified and strengthened the final product. The authors also wish to acknowledge the many individuals at Stantec Consulting Ltd., Golder Associates Ltd. and Blackport Hydrogeology Inc. for cross-section interpretations and hydrogeologic analysis used in developing the numerical model. The Ontario Ministry of Natural Resources and Ontario Ministry of the Environment are acknowledged for funding the model refinement and calibration tasks undertaken as part of the Region of Waterloo Tier Three Water Budget and Local Area Risk Assessment study.