A methodology for identifying ecologically significant groundwater recharge areas

Abstract

A methodology has been developed to delineate recharge areas that support groundwater-dependent ecological features such as wetlands and coldwater streams. Numerical groundwater models are employed together with particle tracking techniques to link recharge areas to specific ecological features. Kernel density estimation methods are used to identify areas that contain the maximum number of particle endpoints per unit area. These high-density endpoint clusters are subsequently mapped as ecologically significant groundwater recharge areas (ESGRAs). The technique can be utilized with either finite-element or finite-difference groundwater models which are loosely coupled or fully integrated with a distributed hydrologic model. A good representation of both the surface water and the shallow groundwater system is required and the models should be based on a rigorous understanding of the regional, watershed and local-scale geologic, hydrogeologic and hydrologic setting. This manuscript describes how the methodology was employed to delineate ESGRAs for the Oro North, Oro South and Hawkestone Creeks subwatersheds, located in the northwest portion of the Lake Simcoe watershed in southern Ontario, Canada. The three study subwatersheds are connected to the Oro Moraine complex, a high-recharge feature that supports multiple watersheds and the regional groundwater flow system. A transient, integrated surface water/groundwater model was developed for the watersheds flanking the moraine and was calibrated to daily streamflow and groundwater-level observations. Particle tracking was then used to define the groundwater flow system between all streams and wetlands and the areas of recharge, and provided insight into the hydrologic function of the Oro Moraine. After optimization of the cluster analysis procedure, 24% of the 125-km² study area was mapped as having recharge areas that support significant ecological features. As required by the Lake Simcoe Protection Plan, the ESGRAs will be protected from future development to ensure the maintenance of the groundwater-fed ecosystems they support. The methodology described in this paper provides a consistent, objective and technically sound means of identifying and delineating ESGRAs, and has recently been applied to other watersheds in southern Ontario.

Une méthode employant des modèles d'écoulements souterrains avec des techniques de traçage de particules a été développée afin de relier des zones de recharge aux écosystèmes qu'elles alimentent. Des méthodes d'estimation de densité sont utilisées pour identifier les zones qui contiennent le nombre maximum de points de terminaison de particules par unité de surface. Ces zones de haute densité sont ensuite cartographiées comme des zones de recharges significatives pour les écosystèmes (ZRSE). Cette technique peut être utilisée avec des modèles d'écoulements souterrains utilisant soit les éléments finis ou les différences finies qui sont couplés ou intégrés avec un modèle hydrologique distribué. Une bonne représentation des écosystèmes et des eaux souterraines peu profondes doit être incluse dans le modèle. Ce manuscrit décrit comment la méthodologie qui a été utilisée pour délimiter les ZRSEs pour les sous-bassins d'Oro Nord, Oro Sud et Hawkestone situés dans la partie nord-ouest du bassin du lac Simcoe, dans le sud de l'Ontario, Canada. Les trois sous-bassins hydrographiques de l'étude sont reliés à la moraine d'Oro, une zone de forte recharge qui alimente plusieurs bassins versants et le système régional d'écoulements souterrains. Un modèle intégré a été développé pour les sous-bassins reliés à la moraine et a été calé de manière à reproduire les débits journaliers en rivière et les niveaux piézométriques. Ensuite les particules ont été utilisées pour définir le système d'écoulement des eaux souterraines entre les écosystèmes et les zones de recharges. Après l'optimisation des paramètres d'analyses en grappe, 24% de la zone d'étude (125 km²) a été définie comme zone de recharge alimentant les écosystèmes. Selon les requis du Plan de protection du lac Simcoe, les ZRSEs seront protégées contre le développement futur afin d’assurer le maintien des écosystèmes reliés aux eaux souterraines. Enfin, la méthodologie présentée ici est cohérente, objective et techniquement fiable pour identifier et délimiter les ZRSE. Elle est présentement appliquée à plusieurs sous-bassins versants du sud de l'Ontario.

Introduction

Lake Simcoe is a 722-km² lake located in southern Ontario, Canada, 100 km north of Toronto. Lake Simcoe’s 2840 km² watershed has been largely converted to agricultural land; however, the watershed’s urban areas are growing, exerting greater pressure on the remaining natural features including provincially significant wetlands, woodlands and coldwater streams. Consequently, human activity in the area has proven detrimental to the ecological function of Lake Simcoe itself due to excessive nutrient loading from agricultural runoff, pollutants and invasive species (Ontario Ministry of the Environment [MOE] 2009). To address this concern, the Lake Simcoe Protection Act was passed by the Ontario Legislature in 2008. Subsequently, in 2009, the MOE released the Lake Simcoe Protection Plan (LSPP) designed to protect, improve and/or restore the ecological integrity of the Lake Simcoe Watershed and its key natural heritage features and functions (MOE 2009).

The objectives of the LSPP relate to the protection, improvement and restoration of hydrologic elements that contribute to the ecological health and function of key natural heritage features within the Lake Simcoe watershed, and aim to restore a self-sustaining coldwater fish community in Lake Simcoe. Several policies within the LSPP aim to define, identify and protect ecologically significant groundwater recharge areas (ESGRAs). Through Policy 6.37-SA, ESGRAs are defined as areas of land that are responsible for supporting hydraulic pathways that sustain sensitive groundwater-dependent ecosystems (GDEs) such as coldwater streams and wetlands. The connection between ESGRAs and significant stream and wetland features is sketched in Figure 1. To establish the ecological significance of a recharge area, the linkages between the recharge area and an ecologically significant feature (e.g. a reach of a stream, a wetland, or an area of natural or scientific interest) must be identified. Once identified, LSPP policy 6.38-DP requires municipalities to incorporate ESGRAs within their official plans with policies to protect, improve or restore the quality and quantity of groundwater in these areas, and the function of the recharge areas.

Figure 1. The role ecologically significant groundwater recharge areas (ESGRAs) play in supporting ecologically significant features (reproduced with permission from Lake Simcoe Region Conservation Authority 2014).

The LSPP guidance documents do not specify a method for establishing the linkage between recharge areas and ecological features required to delineate ESGRAs. The purpose of this paper is to provide a description of the ESGRA methodology developed by Earthfx (2012) for the Lake Simcoe Region Conservation Authority (LSRCA). The ESGRA methodology is described below by use of a modelling study commissioned by the LSRCA (Earthfx 2013a, 2013b). Since its development, the ESGRA methodology described herein has been successfully applied to other southern Ontario watersheds (GENIVAR 2013; Matrix Solutions Inc. 2013b; Earthfx 2014a, 2014c).

Background

Groundwater recharge can be defined as the component of precipitation that percolates through the unsaturated zone and replenishes the groundwater system at the water table, while groundwater discharge is that portion of groundwater removed from the saturated zone by drainage or exfiltration (Freeze and Cherry 1979). Focused groundwater discharge at the surface can be identified as seeps and springs, which are critical to maintaining GDEs in wetlands and streams. According to Batelaan et al. (2003), “the relationship between groundwater recharge and discharge is one of the most important aspects in the protection of ecologically valuable areas (86).

The reliance of wetlands and streams on the groundwater system is quite broad, including the maintenance of aquatic habitat (Imhof et al. 1996; Power et al. 1999; Malcolm et al. 2004; Vrdoljak and Hart 2007), maintaining contiguous buffer zones that promote stream productivity by maintaining baseflow discharge to support fish migration (Beatty et al. 2010), influencing water quality by augmenting dissolved oxygen and providing a thermal refuge for aquatic species (Power et al. 1999; Hayashi and Rosenberry 2002; Malcom et al. 2004), and maintaining vegetative cover in streams (D.M. Merritt and Bateman 2012) and riparian areas (Kuglerová et al. 2014). These positive benefits provided by groundwater discharge features can be impaired either by the alteration of the landscape or through groundwater withdrawals (Curry et al. 2002; Batelaan et al. 2003; Webb and Leake 2006; Waco and Taylor 2010; Menció and Mas-Pla 2010).

Identifying the pathways between GDEs and recharge areas has been recognized as an important field of research, yet little practical work has been undertaken at the watershed scale (Johansen et al. 2014). Much of the work in quantifying groundwater discharge to ecological systems includes the development of tools used to characterize the in-stream flow regime through statistical, geomorphological and biological indicators (Richter et al. 1996; Annear et al. 2004), the use of stream-flow records to identify drought conditions and their relationship to human development and ecological systems (Tallaksen and van Lanen 2004), the quantification of stream-flow components by simplified baseflow separation routines (Eckhardt 2008) or by more complex approaches using integrated surface–subsurface hydrological models (Partington et al. 2013), and through detailed feature-scale field investigation (Gehrels and Mulamoottil 1990; Roulet 1990).

These methods, however, tend to be reactive, serving only to identify alteration of the groundwater discharge regime after some stress has been imposed, while the LSPP policy 6.38–DP specifically seeks a proactive approach by means of identifying lands that require protection from development (i.e. potential stress) to maintain the current rates of recharge to sustain the GDEs they support. Due to the nature of typical groundwater systems, travel times from the point of entry (i.e. recharge) to the point of discharge can be on the order of tens to thousands of years, precluding field methods as a means of establishing these linkages. Therefore, to meet the requirements of the LSPP, the methodology must rely on numerical models, which have shown promise in linking the various hydrological and hydrogeological processes at the corresponding spatiotemporal scales required to identify recharge-to-discharge pathways (Batelaan et al. 2003; Hunt et al. 2013; Johansen et al. 2014).

Methodology

Guidance documents for the LSPP (LSRCA 2014) define ESGRAs as “areas of land that are responsible for replenishing groundwater systems that directly support sensitive areas like coldwater streams and wetlands” (4). To establish the ecological significance of the recharge area, a linkage must be established between the recharge area and the ecologically significant feature. The methodology described below is based on the investigation performed by Earthfx (2013a, 2013b) for the LSRCA.

The delineation of ESGRAs first proceeded by mapping sensitive features of interest identified by the LSRCA, where groundwater discharge is likely occurring over the long term. From this mapping, the regional hydrogeological setting was conceptualized to identify the likely sources of recharge supporting the identified GDEs. Second, a numerical model was constructed, large enough to encompass all potential pathways while maintaining a resolution fine enough to capture the geometry of the ecological features. Lastly, the three-dimensional groundwater flow fields produced by the model were used to determine flow pathways. Based on these pathways, ESGRAs were determined by connecting the points of recharge (ESGRA) to the points of discharge (GDE).

Study area

The Oro North, Oro South and Hawkestone Creeks subwatersheds (Earthfx 2013a) are located in the northwest portion of the Lake Simcoe watershed, as shown in Figure 2. These largely undeveloped subwatersheds drain the southeast flank of the Oro Moraine, an east–west trending feature that rests on the large till upland (Beckers and Frind 2000). The Oro Moraine lies mainly in the Simcoe uplands physiographic region (Chapman and Putnam 1984), which is fringed on the south, west and east by the Simcoe Lowlands, and partly on the north by a segment of the Carden Plain. The Simcoe uplands are dominated by till plains and broad erosional valleys that contain either sand or clay plains. Chapman and Putnam (1984) classify the till plains in the uplands as drumlinized till plains, but there are few drumlins mapped on the uplands in the study area. Lacustrine sand plains predominate in the low-lying areas of the study area, with some clay plain in the northern part of the study area.

Figure 2. Oro Moraine model extents with study subwatersheds.

A number of wetland features identified as significant are present within the study area, especially along the flanks of the Oro Moraine and in the low-lying valleys. Land-use inventory mapping was reviewed to capture these hydrologic features using the Southern Ontario Land Resource Information System (SOLRIS, version 1.2) (Ontario Ministry of Natural Resources [OMNR] 2008). Marshes, bogs, swamps and fens represent 15.5% of the study area. Open water represents 1.1% of the study area, mostly reflecting two small lakes: Bass Lake and Little Lake.

The Oro Moraine has an area of 165 km² and is composed mainly of sand and gravel. It is bounded to the north and west by sub-glacial tunnel channels that were cut through the underlying deposits and later infilled (Burt and Dodge 2011). There are a number of abandoned beaches in the study area, which developed along the shorelines of postglacial lakes, in the low-lying areas and on the flanks of the till uplands (Barnett 1988).

The mapped surficial geology of the area west of Lake Simcoe has been incorporated into the Ontario Geological Survey (OGS) digital compilation map of southern Ontario Quaternary geology (OGS 2010). Drift in the model area is generally quite thick and ranges from 0 to 250 m. In the lowland areas it usually ranges from about 50 to 100 m, with thin (< 15 m thick) drift in the northern part of the model area, particularly to the east between Lake Simcoe and Lake Couchiching. Overburden thickness in the till uplands is typically 60 to 175 m, with the uplands of the Oro Moraine exhibiting a thickness of 250 m (Burt and Dodge 2011).

Land surface topography was based on the 5-m digital elevation model produced by the OMNR. Higher elevations occur along an east–west ridge formed by the Oro Moraine in the center of the study area, with the highest elevation at 405 m above sea level. Local relief ranges from 20 m to more than 150 m on the north side of the moraine. Areas of hummocky topography occur on top of the Oro Moraine and act to prevent surface runoff and focus infiltration (Earthfx 2013b).

Monthly climate normals were derived from daily historical climate data provided by Environment Canada. Monthly average temperature ranges from approximately –8.3°C in January to 20.2°C in July. Temperature data are consistent among the climate stations in the study area. Monthly rainfall rates for the stations are similar, although rainfall exhibits an increasing trend to the east with the total annual rainfall averaging 711 mm/yr. Snowfall rates are more variable among stations, averaging 269 mm/yr snow water equivalent. Total precipitation based on these data averaged 980 mm/yr. Precipitation is higher from August to January, averaging 93 mm/month, and lowest from February to April, averaging 63 mm/month.

Model development

The OGS previously developed a 23-layer hydrostratigraphic model of the moraine and surrounding tunnel channel features (Burt and Dodge 2011). The model featured two bedrock layers and 21 Quaternary (overburden) layers, comprised of alternating aquifer and aquitard units. The OGS geologic surfaces formed the basis of an integrated groundwater–surface water model developed with the United States Geological Survey (USGS)’s GSFLOW code (Markstrom et al. 2008). GSFLOW is based on a linkage between two sub-models: the Precipitation Runoff Modelling System (PRMS) hydrologic sub-model (Leavesley et al. 1983) and the MODFLOW-NWT groundwater sub-model (Niswonger et al. 2011). GSFLOW has been utilized to investigate surface water–groundwater interactions in a number of recent studies (Huntington and Niswonger 2012; Niswonger et al. 2014; Tanvir Hassan et al. 2014; Woolfenden and Nishikawa 2014).

The modelled area was extended significantly beyond the three study subwatersheds to include the entire Oro Moraine as well as several adjacent watersheds (Figure 2). This was done to capture the hydrologic and hydrogeologic function of the moraine as a whole, as well as to ensure that the model was constrained by realistic boundary conditions. This extension had the added benefit of increasing the data available for model calibration. The total modelled area was 805 km² and incorporated 1072 km of streams and 87 lakes and ponds, bounded to the east by Lake Simcoe and Lake Couchiching. The OGS conceptual model was reduced to 12 layers when incorporated into the groundwater sub-model that was constructed at a 100-m resolution (343 columns by 350 rows). As shown in Figure 3, a northwest–southeast section through the study area, the geologic layering is highly irregular resulting in complex three-dimensional groundwater flow patterns. The distributed hydrologic sub-model was constructed on a 50-m grid resolution (322,152 unique hydrologic response units) and parameterized with detailed land use and surficial geology mapping discussed above.

Figure 3. Section A–A’ illustrating the hydrostratigraphic layers incorporated into the GSFLOW model overlain with forward particle tracks from the Oro Moraine toward surrounding ecological features.

The hydrologic sub-model computes canopy interception, snowmelt, infiltration, overland runoff, actual evapotranspiration, interflow and recharge on a daily time step. Climate inputs to the model include hourly NEXRAD precipitation, minimum and maximum air temperature, and global solar radiation. Significant feedback between the groundwater and hydrologic models occurs in areas where the water table is shallow and intersects the soil zone limiting recharge. Interaction also occurs where groundwater discharges to (or is recharged by) lakes, wetlands or streams. All mapped streams, lakes and wetlands were represented explicitly in the GSFLOW integrated model. Streams are represented with the stream-flow routing module SFR2 (Niswonger and Prudic 2005), which routes flow and calculates leakage between streams and the underlying aquifer, and the lake simulation module, LAK3 (M.L. Merritt and Konikow 2000) is used to represent lakes and shallow wetlands.

Observed groundwater levels (also referred to as aquifer heads or potentials) served as the primary calibration targets for the groundwater flow model. The groundwater sub-model was first calibrated in steady state against 3436 static water levels obtained from the Ministry of Environment and Climate Change Water Well Information System (WWIS). Aquifer and aquitard properties were adjusted within reasonable limits during the model calibration to improve the match. The steady-state simulation produced a mean error of 0.3 m with a root mean squared error of 8.1 m when compared to the WWIS data (Figure 4b).

Figure 4. (a) Simulated model and observed streamflow at Water Survey of Canada Coldwater River at Coldwater (02ED007). (b) Scatter plot of simulated versus observed groundwater heads. (c) Simulated and observed heads at PGMN Well #W0000442. WY: water year.

The transient integrated model run at daily time steps was calibrated for water years (WY) 2005 to 2009 inclusive, the period with the greatest amount of observation data. The model was calibrated to match observed daily streamflow and estimated baseflow at four Water Survey of Canada stream gauges, as well as transient groundwater observations. A daily Nash–Sutcliffe efficiency (NSE, Nash and Sutcliffe 1970) factor of 0.63 was achieved at the primary streamflow calibration gauge, Coldwater River at Coldwater (02ED007), shown in Figure 4a. The model calibration matches the WY2007 extreme summer low flows, as well as the higher baseflow and event peaks observed during WY2008 and WY2009. Runoff volumes also compared favorably, with the model achieving a monthly NSE of 0.80. Groundwater heads were further refined in transient mode against observed continuous groundwater levels at eight provincial groundwater monitoring network (PGMN) wells in the model area. The calibration to PGMN Well W0000442 located on the northwest flank of the Moraine is shown in Figure 4c.

Upon satisfactory calibration of the various sub-models and the integrated transient GSFLOW model, a 25-year simulation was performed at a daily time step. The generated long-term average net recharge (recharge less discharge) from the numerical model was employed to delineate the ESGRAs.

Particle tracking

Particle tracking is an accepted methodology for visualizing groundwater flow paths predicted by numerical modelling (Buxton et al. 1991; Zheng and Bennett 1995; Modica et al. 1998; Rock and Kupfersberger 2002; Batelaan et al. 2003; Matula et al. 2014), and is particularly useful in areas with complex, three-dimensional flow fields. Particle tracking can be accomplished with a number of numerical modelling packages. For example, the ESGRA delineation methods described herein were initially tested by Earthfx (2012) with FEFLOW version 6.0 (DHI-WASY 2010).

To initiate tracking, virtual particles are released at points within selected model cells and then tracked in the direction of flow until the particles exit the model domain. For forward tracking in the direction of flow, particles are usually introduced in a uniform distribution across the model domain (Batelaan et al. 2003). The particles are then traced to a point of groundwater discharge at surface or to where they exit through model boundaries. With reverse (or backward) particle tracking, particles are introduced in a dense distribution at the points of known groundwater discharge or within (for example) an ecological feature of interest, and are tracked backward in time.

In Batelaan et al. (2003), a regional groundwater flow system was simulated using the MODFLOW code (Harbaugh and McDonald 1996), and particle tracking techniques were used to delineate the recharge area contributing to each ecologically valuable discharge area. A forward tracking approach was employed by tracking particles from every model cell and delineating the capture zones associated with each individual discharge area. While forward tracking-based capture zone analysis is often used to assess groundwater vulnerability (Snyder et al. 1998), the method is unsuitable for ESGRA delineation. The requirement of LSPP 6.38-DP is to define areas of significant recharge that directly sustain GDEs, rather than determining areas that could potentially contribute to the GDEs.

With reverse particle tracking, specific ecological features can be analyzed in detail as particles released from GDEs will converge on specific areas where the contributing recharge is likely to occur. The density of reversed particle endpoints can then be used as a surrogate indicator of the significance of the recharge areas in terms of their connectivity to GDEs.

Particle tracking analysis for the Oro Moraine model was undertaken with the USGS program MODPATH (Pollock 2012). The groundwater velocity flow field used for the particle tracking analysis was produced by running the Oro Moraine model in steady state. The steady-state model was supplied with a long-term average net recharge simulated by the integrated Oro Moraine model over a 25-year period. Particle tracking analysis performed using steady-state velocity fields overcomes limitations associated with misrepresenting wetland or lake bathymetry (Trigg et. al. 2014).

Environmentally sensitive features including coldwater streams and significant wetlands were identified by the LSRCA and are shown in Figure 2. Using MODPATH, the particles were released along the top face of model cells from within these features. A 5 × 5 m grid of particles was released within wetland boundaries and over stream segments. In order to capture interactions with the stream channels and the riparian areas adjacent to the stream, 400 particles were released in a 100 × 100 m GSFLOW model cell containing a stream segment. Particle path lines were produced by MODPATH reflecting the connection between ecological features and the originating area of recharge that supports the features. A simplified sketch of particle path lines tracked backwards from a study area feature (the Bluffs Creek West wetland located in North Oro) is presented in Figure 5.

Figure 5. Example of particle-tracking from a significant study area feature (Bluffs Creek West wetland, Ontario, Canada) backward to areas of local recharge. ESGRAs: ecologically significant groundwater recharge areas.

Cluster analysis

The reverse particle tracking exercise produces a set of coordinates associated with the location at which the particles tracked in reversed time exited the model domain at the land surface; these coordinates are deemed endpoints. With the exception of certain unlikely situations (such as a perfectly uniform flow field or a flow field that extends perfectly radially with uniform velocity from a discharging feature), reverse-particle path lines will tend to converge, resulting in endpoints that have a tendency to cluster. ESGRA analysis is thus based on the assumption that clusters of greatest endpoint density are an indication of areas on the landscape whose recharge has a greater propensity of contributing directly to the GDEs from which the particles were released.

A simple cluster analysis approach that was first investigated by Earthfx (2012) involved overlaying a grid of a given resolution and counting the number of endpoints contained within every grid cell. Although this approach was able to provide a rough estimate of potential ESGRAs, the resulting ESGRA coverage was ultimately dependent on grid spacing, origin and orientation. Furthermore, an endpoint that lies near the edge of a grid cell would only have influence within the particular cell it falls in, and not its spatial location relative to neighbouring grid cells. A solution for delineating ESGRAs that avoided these limitations was accomplished by performing the cluster analysis using a bivariate kernel density cluster analysis technique to quantify endpoint density.

Details of the method developed to objectively evaluate endpoint clusters and delineate ESGRAs are presented in the Earthfx (2012) study. The method proposed was adopted from published, peer-reviewed cluster analysis methodologies. This technique was further tested and refined so that it could be applied to other sub-watersheds, ensuring that the delineation of ESGRAs can be conducted in a consistent manner. The cluster methodology is based on a non-parametric bivariate symmetric kernel density estimation technique modified from Wand and Jones (1993) as:(1)

where is the Euclidian distance between an observation point and endpoint , n is the total number of endpoints and h is the kernel bandwidth. In the context of ESGRA cluster analysis, by first selecting an appropriate kernel bandwidth (h), the contribution that a selected endpoint () makes to the relative density at observation point () can be computed. Once all endpoints (i = 1, …, n) have been evaluated, the sum of their contributions represents the expected endpoint density () at the observation point.

A final pass of the bivariate kernel density cluster estimation field () is applied to eliminate ESGRAs of small areal extent, infill holes present in an ESGRA contours and eliminate areas of relatively small density. In Earthfx (2012), the minimum allowable ESGRA was set to 0.045 km², which corresponded to the average model element area. Similarly, holes in the ESGRA contours of less than 0.045 km² were filled in to produce continuous ESGRA delineations. Areas of small densities were removed by first defining a delineation cut-off threshold () and eliminating all areas where the calculated density is less than a ^th of the maximum evaluated (i.e. eliminate all areas where ). Typically the delineation cut-off threshold is on the order of , meaning that any area where the evaluated is less than 1% of the maximum evaluated density () is removed from the final ESGRA coverage.

Results and discussion

A total of 1.6 million particles were released from select areas of ecological significance within the model. Particle tracking via MODPATH resulted in a data set of reverse particle tracking endpoint coordinates. Endpoint cluster analysis was optimized by permuting the bandwidth (h) and cut-off threshold () of the bivariate kernel density estimator. Optimization was achieved by determining the greatest number of endpoints captured by the smallest possible area. Table 1 presents the number of endpoints contained within a cluster proportional to the total number of particles released. Table 2 presents the corresponding percent of the total study area delineated as a cluster for a given combination of kernel density estimator parameters.

Table 1. Percent of endpoints covered by potential ecologically significant groundwater recharge areas (ESGRAs) with varying smoothing parameter () and delineation threshold ().

ε	h = 10 m	h = 25 m	h = 50 m	h = 100 m
10	7.3%	15%	29%	54%
20	18%	37%	61%	61%
100	65%	88%	97%	97%
200	83%	96%	99%	99%
1000	99%	100%	100%	100%

Table 2. Percent of study area covered by potential ecologically significant groundwater recharge areas (ESGRAs) with varying smoothing parameter () and delineation threshold ().

ε	h = 10 m	h = 25 m	h = 50 m	h = 100 m
10	0.08%	0.39%	1.8%	8.4%
20	0.36%	2.0%	7.4%	23%
100	5.3%	16%	31%	55%
200	11%	24%	40%	63%
1000	24%	37%	53%	75%

Of the particles released, approximately 955,200 (60%) were released into discharging cells. These were used for endpoint analysis and ESGRA delineation. The remaining particles were released into cells that were found to be locally recharging the groundwater system (e.g. a losing stream reach or portions of an identified wetland disconnected from the phreatic surface). These particles did not leave their starting cells and were therefore excluded from the endpoint analysis. Of the valid endpoints, about 894,700 (93.7%) remained within the study area subwatersheds (i.e. Oro North, Oro South and Hawkestone Creeks subwatersheds), while the remaining particles tracked beyond subwatershed boundaries and terminated in areas further west along the Oro Moraine. While the number of particles exiting the study subwatersheds was not large, it did indicate that some features, in particular the headwaters of Hawkestone Creek, are likely receiving lateral groundwater inflow from outside the subwatershed.

Based on these tables and from visual inspection of the resulting cluster maps (such as Figure 5), the kernel density bandwidth (h) was set to 25 m and a delineation threshold to . With these parameter values, the clusters were deemed ESGRAs because they consistently met the following criteria:

(1)	Rejection of endpoints that clearly did not belong to any cluster;
(2)	Delineation of clusters with a relatively high density of particle endpoints; while
(3)	Avoiding spreading into areas where endpoints did not exist.

The particle endpoint distribution was mapped onto a 25 × 25 m uniform gridded surface to produce the final ESGRA mapping (Figure 6). Of the particles released from discharging cells, 96% are included in ESGRA zones.

Figure 6. Delineated ecologically significant groundwater recharge areas (ESGRAs) contributing to significant ecological features within the study area.

The subwatershed area covered by ESGRAs varied among the study subwatersheds, with 15% coverage within the Oro South subwatershed versus 26% within Hawkestone Creek subwatershed and 23% within the Oro North subwatershed. The smaller percentage of ESGRA coverage in the Oro South subwatershed is a result of the silty to sandy silt Newmarket till at the surface, reducing recharge relative to the Hawkestone and Oro South subwatersheds.

The particle tracking exercise demonstrated the connectivity between GDEs and specific recharge features (the delineated ESGRAs). Some GDEs were connected to large-scale regional features while others were connected to smaller, local features. For example, the headwaters of Hawkestone Creek exhibited a direct connection to the Oro Moraine, an area of regionally high recharge and significant groundwater storage, whereas along the lower reaches of Hawkestone Creek, GDEs were connected areas of local recharge through adjacent thin sand deposits on till (such as the West Bluffs wetland shown on Figure 5). Simulations were undertaken with the transient integrated model (Earthfx 2013b) to investigate the role of groundwater storage on a watershed scale. Results with the transient GSFLOW model showed that GDEs connected to areas of high groundwater storage proved to be resilient to simulated historical drought conditions, while features fed by local recharge were found to be extremely sensitive.

Forward particle tracking can be further employed to investigate the movement of groundwater from specific features at various spatial scales. Figure 3 shows path lines on a section from particles released over the Oro Moraine and forward tracked to areas of discharge. As can be seen, flow is mainly vertical through the aquitards and lateral through the aquifers. Groundwater moves radially outward from the Oro Moraine and primarily discharges to headwater streams and wetlands along the moraine perimeter. A subset of the flow paths penetrate through the till confining units and flow laterally through the deeper aquifers. The path lines illustrate the direct linkage between the Oro Moraine and the Hawkstone Creek wetlands.

The method was demonstrated using a calibrated GSFLOW model built on a detailed analysis of the geology, hydrostratigraphy, land use and surface topography. The model also incorporated a detailed representation of surface water features (e.g. headwater and first-order streams, wetlands and water bodies) and included neighbouring sub-watersheds to exclude model boundary effects. The method could have been employed with a standalone steady-state groundwater model, but this approach would have required an estimation of distributed recharge. By applying an integrated surface water–groundwater model, uncertainty associated with the estimation of distributed recharge is reduced as groundwater feedback has been accounted for.

The application of a steady-state model does have its limitation in that the flow fields produced are not representative of actual groundwater flow paths. The development of a set of flow fields from a transient model run could provide insight into the seasonal variation of the hydrogeologic system, the reduced/increased contribution to GDEs under drier/wetter conditions and the overall response of the watershed to development and climate change.

Conclusion

A description of a method for producing a quantitative and repeatable technique for the identification of ecologically significant groundwater recharge areas (ESGRAs) was presented. The goal of ESGRAs is to identify regions on the landscape that potentially have a direct hydraulic connection to groundwater dependent ecosystems. The purpose of this methodology is to map areas where existing hydrological function must be maintained in order to sustain the ecological function of systems dependent on groundwater discharge.

The methodology employs reverse particle tracking to link significant ecological surface water features to contributing groundwater recharge areas, and cluster analysis to identify concentrations of particle endpoints. Recharge significance is assumed to be proportional to endpoint density and the cluster analysis is optimized to preferentially contour areas of the highest endpoint density. The resulting ESGRA map can be used to visualize and compare these areas on a watershed-by-watershed basis, under both existing or future development, and hypothetical climatic conditions. Reverse particle tracking from known ecological features is a useful and direct means of identifying contributing recharge locations. The technique can be used to make comparisons among adjacent watersheds and ecological features, and between current and proposed land development.

This methodology has direct application to many jurisdictions within Ontario, specifically given the 10-year investment in regionally distributed surface water and groundwater numerical models that was made under the Clean Water Act. These models are highly sophisticated and built to a large spatial scale at a fine resolution, with some incorporating the full integration of atmosphere–surface–subsurface processes (AquaResource 2011; Earthfx 2013a, 2014b; Matrix 2013a, 2014; Golder and Associates Ltd. 2014). With this unprecedented wealth of groundwater management tools in the hands of conservation authorities and municipalities, it is now acknowledged that there is a need to leverage this investment by incorporating them into a groundwater knowledge management system (Holysh and Gerber 2014). The ESGRA methodology outlined in this paper is transferable to many models created under Source Water Protection, and could be used to gain additional value from these existing modelling efforts.