Using water stable isotopes for tracing surface and groundwater flow systems in the Barlow-Ojibway Clay Belt, Quebec, Canada

Abstract

This study aims to improve the understanding of surface and groundwater flow systems based on water stable isotope data in a 19,549 km² region of the Barlow-Ojibway Clay Belt, in western Quebec, Canada. The available geochemical database contains 645 samples including precipitation, snow cores, surface waters, groundwater and springs. All samples were analyzed for water stable isotopes (δ²H-δ¹⁸O) and complementary tritium analyses were conducted on 98 groundwater and spring samples. Precipitations depict a clear temperature-dependent seasonal pattern and define a local meteoric water line (LMWL) without a latitudinal trend in δ²H-δ¹⁸O. Samples collected from the snowpack plot on the LMWL, suggesting that the bulk snowpack preserves the isotopic composition of precipitation throughout the frozen period, prior to the spring snowmelt. Surface water samples define a local evaporation line (LEL), and evaporation over inflow (E/I) ratios range between 0 and 36%. Groundwater and spring samples are evenly distributed around the LMWL, suggesting that evaporation processes are limited prior to infiltration and that surface waters do not significantly contribute to groundwater recharge. Shallow unconfined aquifers present a greater variability in δ²H-δ¹⁸O compared to confined aquifers located farther down gradient, suggesting the mixing of varied recharge waters along the regional groundwater flow system. A three-component mixing model based on isotopic and specific electrical conductivity data allows the quantification of such mixing processes. The interpretation of isotopic data constrains a regional-scale conceptual model of groundwater flow systems and describes processes related to the timing of recharge, evaporation, mixing and discharge.

Cette étude vise à mieux comprendre les systèmes d’écoulement dans une région de 19 549 km² au sein de la ceinture argileuse Barlow-Ojibway (ouest du Québec, Canada). La base de données contient 645 échantillons analysés pour les isotopes stables de l’eau (δ²H-δ¹⁸O) et inclue les précipitations, la neige au sol, les eaux de surface, les eaux souterraines et les sources alors que 98 échantillons d’eau souterraine et de sources ont été analysés pour l’activité tritium. La composition isotopique des précipitations présente un cycle saisonnier dépendant de la température et définissant une droite locale des eaux météoriques (DLEM) sans tendance δ²H-δ¹⁸O reliée à la latitude. Les données de neige au sol tombent sur la DLEM, ce qui suggère que le couvert nival préserve sa composition isotopique en hiver, avant la fonte des neiges printanières. Les échantillons d’eau de surface définissent une droite évaporatoire locale (DEL) et les rapports de l’évaporation sur les flux intrants variant entre 0 et 36 % à l’échelle régionale. Les données d’eau souterraine et de sources suggèrent que les processus d’évaporation sont limités avant l’infiltration et que les eaux de surface contribuent peu à la recharge. Les aquifères à nappe libre peu profonds présentent une plus grande variabilité isotopique par rapport aux aquifères à nappe captive situés en aval le long des lignes d’écoulement, ce qui suggère une atténuation de la variabilité isotopique issue de la recharge le long des systèmes régionaux. Un bilan géochimique appuyé sur les données isotopiques et de conductivité électrique de l’eau permet de quantifier ces mélanges. L’interprétation des données isotopiques permet de contraindre des modèles conceptuels représentant les systèmes d’écoulement régionaux et de documenter les processus de recharge, d’évaporation, de mélange et de résurgence.

Introduction

Groundwater plays key functions in several ecosystems and represents an essential source of potable water for approximately half of the world’s population (Taylor et al. 2013). Nevertheless, groundwater resources have been heavily used for human water supply and agriculture for many years on a global scale (WWAP 2009), and climate change will most likely affect the hydrological cycle and aquifers in the near future (Green et al. 2011). Groundwater resources are expected to be increasingly stressed in response to demographic growth and growing water demand (Vörösmarty et al. 2000). It is estimated that by 2025, nearly 50% of the world’s population will be living in water-stressed regions (WHO 2017). Such considerations stress the critical need to better understand hydrogeological processes at various scales to support concrete actions aimed at protecting groundwater resources, to maintain ecosystems and meet water supply needs of future generations. Regionally and locally, groundwater protection must rely on land management strategies that are based on an appropriate understanding of the hydrogeological system. This is especially important in areas where most of the population relies on groundwater for potable water supply, as is the case for the Abitibi-Témiscamingue region in western Quebec. This region is renowned for its large eskers and moraines and for the remarkable quality of groundwater found within the aquifers associated with these formations (Veillette et al. 2004; Collini et al. 2007). However, these shallow unconfined aquifers are highly vulnerable to contamination and face increasing human pressures such as forest operations, sand and gravel extraction, former in-trench disposal sites, private and municipal wells and commercial water bottling (MRNF 2006; Nadeau et al. 2015; Cloutier et al. 2016). The impacts on the groundwater resources associated with eskers might also affect aquifers located down-gradient, further highlighting the need to better understand regional hydrogeological processes. In this context, providing regional stakeholders with a coupled understanding of physical and geochemical conditions is mandatory to ensure that land-use decisions related to groundwater protection are based on reliable information.

Geochemical tracers are included in a wide spectrum of hydrological studies. They are increasingly used to delineate flowpaths and quantify the processes affecting surface and groundwater resources at various scales and in various regions of the world (e.g. Kortelainen and Karhu 2004; Wassenaar et al. 2011; Jeelani et al. 2013). Among the wide range of geochemical tracers that have been exploited, water stable isotope (²H-¹⁸O) data have provided unique insights into hydrological processes (Araguás-Araguás et al. 2000; Gibson et al. 2005; Birks and Gibson 2009). The stable isotopes of water are partitioned within a hydrological system, mainly in response to systematic fractionation mechanisms occurring within the water cycle, including phase-change (solid–liquid–vapour) processes, diffusion and rainout (Craig 1961; Rozanski et al. 1993; Gat 1996; Clark and Fritz 1997). The isotopic composition of water is widely used in theoretical and applied hydrological studies to document sources of atmospheric moisture and precipitation patterns (Dansgaard 1964; Fritz et al. 1987), snow and ice-melting (Siegel and Mandle 1984), evaporation processes in cold regions (Gat et al. 1994) and drylands (Kebede et al. 2009), sources mixing (Ferguson et al. 2007; Wolfe et al. 2007) and groundwater recharge (Praamsma et al. 2009), among other applications.

Cloutier et al. (2013, 2015, 2016) provided a comprehensive assessment of groundwater resources in the Abitibi-Témiscamingue region, within the framework of the Programme d’Acquisition de Connaissances sur les Eaux Souterraines (PACES) supported by the Quebec Ministry of the Environment. A regional-scale conceptual hydrogeological model representing regional flowpaths from unconfined aquifers, located in recharge areas, to confined aquifers, found beneath the clay belt, was proposed through the course of these studies. Nevertheless, the proposed conceptual model is mainly constructed on the basis of physical hydrogeological data, whereas geochemical data have not been yet fully exploited, which leads to important limitations in the understanding of regional hydrogeological dynamics. Fitting into the pre-established framework, this study aims to improve the understanding of regional surface and groundwater flow systems in the southern portion of the Barlow-Ojibway Clay Belt (Abitibi-Témiscamingue, Quebec, Canada) on the basis of geochemical data. The present study focuses on the stable (²H-¹⁸O) and radioactive (³H) isotopes of the water molecule. The proposed approach involves a regional-scale comprehensive assessment of isotopic data within all hydrological components: precipitations, winter snowpack, surface waters, groundwater and springs. The specific objectives of the present study are (1) to document the atmospheric processes driving spatiotemporal variations in precipitation isotopic composition, (2) to assess the potential interactions and mixing processes between surface waters and aquifers, and (3) to decipher the mechanisms controlling the water isotopic variability in regional aquifers. Ultimately, the data are interpreted to support an improved conceptual model of regional groundwater flow systems representative of observed hydrogeological and geochemical conditions.

Study area

The study area (Figure 1) covers 19,549 km² and corresponds to the municipalized territory of the Abitibi-Témiscamingue region in Western Quebec, Canada. It encompasses the City of Rouyn-Noranda and the Regional County Municipalities (RCM) of La Vallée-de-l’Or, Abitibi, Abitibi-Ouest and Témiscamingue, for a total population of approximately 148,000 inhabitants (MAMROT 2016). It is estimated that more than 70% of this population relies on groundwater as a source of domestic drinking water. The characteristics of the study area are briefly summarized below. Further details about the hydrogeological framework can be found in Nadeau et al. (2015, 2017) and Cloutier et al. (2016), and within the geological maps of the Geological Survey of Canada (Veillette 1986, 1987a, 1987b; Paradis 2005, 2007; Thibaudeau and Veillette 2005).

Figure 1. Study area with the location of sampling sites. The number (n) of available samples is shown brackets in the legend. The limits of the study area correspond to those of previous studies from Cloutier et al. (2013, 2015). The continental water divide between the St. Lawrence and the James Bay basins is also shown.

Physiography and hydrography

The landscape of the study area is relatively flat, with elevations generally ranging between 280 and 320 m. The main positive reliefs correspond to heterogeneously distributed bedrock outcrops, while the main depressions are associated with large lakes (Abitibi, Témiscamingue and Simard; Figure 1). Despite its featureless landscape, the region is crossed by the continental water divide separating surface waters of the St. Lawrence River Basin (southern sector of the study area, Ottawa River watershed) and those of the James Bay Basin (northern sector of the study area, Moose, Harricana and Nottaway rivers watersheds).

Climate and air masses

The study region is characterized by a cold and humid continental climate with damp summers, rather cold and long winters and a marked seasonality. The average monthly air temperature remains below the freezing point from November to March (based on Quebec Ministry of the Environment climate data; 1981–2010). The region is characterized by a significant latitudinal thermal gradient, and a difference of 2.1°C between the average annual air temperature in Témiscamingue (south) and Abitibi (north) was reported by Asselin (1995). This difference is likely imparted by latitude and by the influence of the continental water divide that creates a discontinuity between the northern and southern sectors of the study area. The presence of numerous lakes in a territory characterized by a gentle relief likely favors the development of buffered microclimates. Such a phenomenon is documented in the vicinity of lakes Abitibi and Témiscamingue (Asselin 1995). Cloutier et al. (2015) provided daily air temperatures and vertical inflows based on interpolated data provided by Centre d’Expertise Hydrique du Québec (CEHQ) for the period ranging from 1900 to 2010 (Poirier et al. 2012). Vertical inflows correspond to the fluxes of water from liquid precipitation and from the melting of snow, and therefore to the fluxes of water that are available for infiltration or surface runoff. If the temperature is below the freezing point, vertical inflows can be null even if precipitation events are observed. In the absence of snow on the ground and when precipitation falls as liquid water, vertical inflows correspond to precipitation. Yearly average vertical inflows range between 790.5 and 856.5 mm at the regional scale. The vertical inflows associated with snow falling during the frozen period will mainly be observed during spring snowmelt. The greatest vertical inflows are recorded in April for most of the study area, except for the northernmost sector where maximum values are observed in May, probably due to later thaw (Cloutier et al. 2015). Minimum vertical inflows are recorded in January over the entire region due to temperatures that generally remain below the freezing point. Regional-scale monthly averages of daily maximum and minimum temperatures reveal that the warmest month is July (with daily average minimum and maximum temperatures of 11.2°C and 23.7°C, respectively) and the coldest month is January (with daily average minimum and maximum temperatures of −23.3°C and −10.3°C, respectively). The mean annual precipitation ranges from 875 mm for La Vallée-de-l’Or RCM to 989 mm in the city of Rouyn-Noranda (Quebec Ministry of the Environment climate data).

The main air masses affecting western Quebec are (1) the cold and dry Arctic; (2) the moist and warm Maritime-tropical air that originates from the Caribbean, subtropical Atlantic and Gulf of Mexico, and migrates northwards following the Mississippi and Missouri valleys before being deflected eastwards across the Great Lakes basin; and (3) the northern Pacific westerlies (Bryson 1966; Bryson and Hare 1974; Hare and Hay 1974). The Pacific westerlies are dominant in the continental climate for much of the year (Hare and Hay 1974). Summers in the study area are warm and humid, reflecting the dominance of Maritime-tropical air punctuated by intrusion of mild and dry westerlies. Winters are influenced by Artic air masses that lead to cold and dry conditions (Bryson and Hare 1974). The Abitibi-Témiscamingue region is characterized by weather variability and instability reflecting the contribution of the different air masses. The dominant winds generally originate from the northwest during winter and the southwest during summer. The polar jet stream, corresponding to the limit between the Arctic and Maritime-tropical air masses, crosses the study area along a sinuous west–east pattern at latitudes roughly corresponding to the position of the continental water divide.

Hydrogeological framework

The regional geological framework, as conceptually represented in Figure 2, is described in detail in numerous previous studies (Veillette 1986, 1987a, 1987b; Veillette et al. 2004; Paradis 2005, 2007; Thibaudeau and Veillette 2005; Nadeau et al. 2015, 2017; Cloutier et al. 2016) and it is thus only briefly presented here. Table 1 provides a summary of the hydrological domains (atmosphere, surface, subsurface), components (precipitation [snow, rain], surface water, groundwater, springs) and aquifers (granular [unconfined and confined], fractured rock [unconfined and confined]) that are discussed in this study, in relation with the main geological units of the region. The bedrock of the study area (geological unit A in Figure 2) is composed of a wide variety of igneous, metamorphic and metasedimentary rocks of the Abitibi and Pontiac sub-provinces, both components of the Superior Geological Province. Groundwater flow within the bedrock is controlled by the architecture of structural discontinuities. Cloutier et al. (2016) estimated hydraulic conductivities ranging between 10⁻⁹ and 10⁻⁴ m/s for this unit. However, there are currently insufficient data to propose a regional-scale interpretation of bedrock hydrogeological characteristics based on structural interpretations. Nevertheless, Rouleau et al. (1999) proposed that groundwater flow within the bedrock of the study area is most likely restricted to depths shallower than 75 m where sub-horizontal fractures are sufficiently interconnected to allow significant groundwater flow. Most of the domestic wells of the region are withdrawing groundwater from this unit using wells that are cased through surficial sediments and open across the bedrock.

Figure 2. Regional hydrogeological conceptual model and water sampling scheme used to characterize the isotopic signature of all hydrological system components. (A) Fractured bedrock; (B) till; (C1) eskers and moraines; (C2) glaciofluvial sediments found between eskers and moraines, under the clay plain; (D) fine-grained deep-water glaciolacustrine sediments; (E) sublittoral sands and beach or eolian sediments; (F) organic deposits. The numbers (1 to 8) correspond to the sampled hydrological components and hydrogeological units, as identified in Table 1.

Table 1. Summarized characteristics of the regional hydrogeological framework. The numbers (1 to 8) of the sampled hydrological components and hydrogeological units correspond to those shown in Figures 2 and 8. The till (B) is not an aquitard, nor is it an aquifer in itself, although it can contribute to the flow of overlying or underlying aquifers.

Domains	Hydrological components	Geological units	Hydrogeological units and aquifers
Atmosphere	(1) Precipitation (rain and snow)	–	–
Surface	(2) Snowpack	–	–
	(3) Surface water
	(4) Spring
Subsurface (groundwater)		C1 Eskers and moraines	(5) Unconfined granular aquifer
		E Surface sediments near eskers
		D Glaciolacustrine clay and silt	Regional aquitard
		C2 Sediments underlying clay	(6) Confined granular aquifer
		B Till	–
		A Fractured rock	(7) Unconfined or (8) confined fractured rock aquifer

The oldest unconsolidated unit of the region corresponds to the till associated with the last glaciation (geological unit B in Figure 2). It is recognized as a compact heterometric geological unit characterized by a matrix composed of sand, silt and clay (Cloutier et al. 2015). The hydrogeological characteristics of this unit are poorly documented because it is rarely exploited for groundwater supply, although it is not considered a regional aquitard (Cloutier et al. 2013, 2015). This till unit is overlain in places by glaciofluvial sediments. The coarse-grained glaciofluvial sediments (eskers and moraines; geological unit C1 in Figure 2) of the region were mainly deposited at the emergence of subglacial meltwater within proglacial Lake Barlow-Ojibway (Nadeau et al. 2015). These formations are characterized by relatively high hydraulic conductivities, 10⁻⁶ to 10⁻¹ m/s according to Cloutier et al. (2016), and constitute the most productive unconfined aquifers of the region (Nadeau et al. 2015). Other glaciofluvial sediments, mainly consisting in finer sand and gravel, are also widespread between the main eskers and moraines of the region (geological unit C2 in Figure 2). This unit is generally characterized by hydraulic conductivities that are slightly lower than those of the eskers and moraines (Cloutier et al. 2013, 2015, 2016). The glaciofluvial sediments are covered in places by fine-grained deep-water glaciolacustrine sediments deposited within proglacial Lake Barlow-Ojibway (identified as the Clay Belt in Figure 1; geological unit D in Figure 2). Owing to its structure consisting in centimetre-scale alternating layers of clay and silt (varves), this unit is considered a regional aquitard (Cloutier et al. 2015), although silt horizons may allow preferential sub-horizontal groundwater flow. Sublittoral sands and beach or eolian sediments (geological unit E in Figure 2) overlying the fine-grained deep-water glaciolacustrine sediments are spatially associated with the glaciofluvial deposits, whereas organic deposits (geological unit F in Figure 2) are widespread across the region, with large peatlands often developing on the flanks of eskers and moraines.

Given this geological framework, and based on previous work from Cloutier et al. (2016) and Nadeau et al. (2017), four main types of aquifers are distinguished within the study region (Figure 2) and used to regroup the isotopic data:

(1)	Unconfined granular aquifers, mainly associated with glaciofluvial formations (eskers and moraines; number 5 in Figure 2). These formations are characterized by relatively high hydraulic conductivities and represent the most productive aquifers of the region, providing a high-quality resource to numerous municipal and private wells.
(2)	Confined granular aquifers, in sectors where till and glaciofluvial sediments are covered by the fine-grained deep-water glaciolacustrine sediments (number 6 in Figure 2);
(3)	Unconfined fractured rock aquifers, essentially restricted to sectors characterized by bedrock and till outcrops (number 7 in Figure 2);
(4)	Confined fractured rock aquifers, in sectors covered by fine-grained deep-water glaciolacustrine sediments (number 8 in Figure 2). Most residences in rural sectors of the region withdraw water from private wells installed in this aquifer.

Cloutier et al. (2013, 2015, 2016) provided estimates of the average groundwater recharge rates (GRR) by calculating the water balance for each parcel of the study region on a monthly basis. The calculations were provided on a 100 m × 100 m grid over the study region, using the following equation:(1)

These water balance calculations are based on the vertical inflows (VI) provided by Poirier et al. (2012) and on hypotheses associated with the regional hydrogeological framework for estimating runoff (R), water storage within the unsaturated soil (ΔW_S) and evapotranspiration (ET). Runoff coefficients were estimated according to land use, soil characteristics and surface slopes, among other things. Potential evapotranspiration rates were quantified using the Thornthwaite (1944) equation with climate data from Poirier et al. (2012). The water storage within the unsaturated zone and the available water for transpiration were evaluated for the different soils based on estimates of water retention capacities, wilting points and thickness of the root zone. The calculated recharge rates range between less than 91 mm/yr in sectors where fine-grained deep-water glaciolacustrine sediments are found and more than 300 mm/yr, mainly in sectors corresponding to eskers/moraines and other unconfined granular aquifers. Overall, maximum recharge rates are observed during the snowmelt period whereas the lowest rates are recorded throughout the winter period, when temperatures are below the freezing point.

Methods

Sampling procedures and in situ measurements

Overall, 645 water samples were analyzed for their isotopic composition (δ²H-δ¹⁸O) as part of this study, over a period ranging from 2006 to 2015 (Figure 1; Table 2). Table 2 provides a summary of the available data and associated sampling periods. From this data set, 98 water samples (only springs and groundwater) were also analyzed for their tritium (³H) content, over a period ranging from 2006 to 2011. The procedures associated with the sampling of precipitation, the snowpack, surface waters, springs and groundwater are presented separately below.

Table 2. Summary of available water stable isotope data. The number of available samples and associated timeline are shown. As reported, most samples were collected during the summer period. Further details about the surface water samples are provided in Table 3.

Type of samples	Number of samples (n)	Sampling timeline
Precipitation	164	Monthly samples at four stations 2009–2015
Snowpack	31	March 2014 (n = 31)
Springs	76	Summer 2009 (n = 30)
		Summer–autumn 2010 (n = 41)
		Summer–autumn 2013 (n = 5)
Kettle lakes	14	Summer 2009 (n = 14)
Surface waters	85	Summer-–autumn 2013 (n = 85)
Confined fractured rock aquifer	105	Summer–autumn 2006 (n = 14)
		Summer–autumn 2010 (n = 26)
		Summer 2013 (n = 65)
Unconfined fractured rock aquifers	38	Autumn 2006 (n = 2)
		Summer 2010 (n = 7)
		Summer–autumn 2013 (n = 29)
Confined granular aquifers	53	Summer–autumn 2006 (n = 2)
		Summer–autumn 2010 (n = 15)
		Summer 2013 (n = 36)
Unconfined granular aquifers	79	Summer–autumn 2006 (n = 2)
		Summer 2009 (n = 14)
		Summer–autumn 2010 (n = 25)
		Autumn 2011 (n = 15)
		Summer–autumn 2013 (n = 23)
Total	645	2006–2015

Precipitation monitoring

Monthly composite precipitation samples were collected at four monitoring stations within the study area (Figure 1; Station 1 = Béarn; Station 2 = Montbeillard; Station 3 = Sainte-Hélène-de-Mancebourg; Station 4 = Amos). Samples were collected monthly from July 2013 to 2015 at stations 1-2-3 (Figure 1), three sites included in the Quebec Ministry of the Environment Climate Monitoring Network. Rainfall samples were collected using 10.1-cm inside diameter All Weather P-2000 standard rain gauges from Productive Alternatives. A thin layer (approximately 1 cm) of paraffin oil (Mineral Oil, Light [NF/FCC], Fisher Chemical) was added to the rain gauges to prevent evaporation, and plastic screens were installed at the base of the rain gauge funnels to block debris and insects. Snowfall samples were collected owing to the participation of local operators/observers associated with the monitoring network. These observers are mandated to make measurements of snow (calculated in terms of equivalent water) intercepted by shielded Nipher snow gauges. Snow samples are collected daily at stations 1 and 3 (but not at station 2 owing to logistical constraints), left to melt overnight and subsequently stored in 1-L high-density polyethylene (HDPE) bottles. Daily samples are mixed in order to obtain composite monthly samples. Paraffin oil (Mineral Oil, Light [NF/FCC], Fisher Chemical) is added to the bottles in order to prevent evaporation. The Amos monitoring station (Station 4), in operation from September 2009 to 2015, is located on the roof of the university campus and is not associated with the Ministère du Développement durable, de l’Environnement et de la Lutte contre les changements climatiques (MDDELCC) climate network. At this station, rain and snow samples are collected daily using a 30-cm inner diameter bucket. Snow samples were left to melt overnight at 4°C and subsequently stored in HDPE bottles. Daily samples were combined to produce a monthly composite sample. In all cases, the volume of monthly composite samples was first evaluated in the laboratory, and smaller aliquots were subsequently stored in 60-mL HDPE bottles and kept at 4°C until stable isotope (δ²H-δ¹⁸O) analyses.

Samples from the snowpack

A sampling campaign was carried out in March 2014 in order to retrieve cores from the snowpack accumulated throughout the winter. Sampling sites were selected in the western part of the study region in an effort to collect data in the vicinity of the grid points used by the CEHQ for calculating vertical inflows (Poirier et al. 2012). At each site, a straight trench varying in length from 5 to 6 m was first dug in a direction oriented perpendicular to the wind. Snow cores (n ≥ 5) were subsequently sampled at ≤ 1-m intervals, on the side of the trench facing the dominant wind. Snow cores were collected using 4.2 cm inner diameter × 1.6 m length HDPE tubes allowing the collection of the snowpack over its entire height. Samples were weighed on site using a Mettler Toledo digital scale with a precision of 0.0001 kg in order to calculate equivalent water contents. For each site, the collected snow cores were transferred and combined within wide-neck HDPE bottles in order to obtain > 800 g composite samples. Samples were kept frozen within the tightly sealed wide-neck bottles during the field campaign and subsequently left to melt overnight at room temperature. The melted samples were filtered using 0.45-μm Waterra FHT cartridges to remove soil particles and/or windblown debris embedded in the snow. The filtered samples were transferred into 30–60-mL HDPE containers and stored at 4°C until stable isotope (δ²H-δ¹⁸O) analyses.

Surface waters

Surface water samples were collected in July 2009 and from July to October 2013. Rivers and streams were sampled near the central part of the main channel, at depths of 0.5–1.5 m below surface. In situ parameters (pH, dissolved oxygen, conductivity, redox potential and temperature) were measured using a YSI 556 MPS probe, a few centimetres downstream of the sampling point. A 7.32-m NASCO Swing Sampler telescopic pole equipped with a 1-L HDPE bottle was used to collect samples from the banks of small rivers and streams. The perch and its container were rinsed on site at a point located a few metres downstream of the sampling site. Whenever sampling from a bridge was possible, the samples were collected using 1-L Waterra ecobailers. A nylon wire equipped with a float was attached to the bailer to indicate the sampling depth. The bailers were systematically rinsed on site, with water from the sampled rivers, prior to sampling. Some of the streams and rivers were sampled immediately at the outlet of large lakes (Table 3). These samples are used to provide an instantaneous estimate of the isotopic composition of water in the epilimnon at the outlet of the lake. Grab samples were also collected from 16 additional lakes, including 14 kettle lakes (Table 3). The samples from Lake Témiscamingue were collected from a boat, at a depth of 1 m using a NASCO Swing Sampler telescopic pole equipped with a 1-L HDPE bottle. The samples from Lake Tee were collected from a boat, approximately 1 m below the surface, using a plastic tube equipped with an inertial valve. The samples from kettle lakes were collected from the shoreline, approximately 1 m below surface, using a NASCO Swing Sampler telescopic pole equipped with a 1-L HDPE bottle.

Table 3. Summary of surface water sampling conditions. Most surface water samples were collected at a shallow depth below the surface, mainly during the summer period. Further details about the sampling procedures are provided in the text.

Surface waters	Names	Sampling techniques	Dates
Kettle lakes	–	1 m from the lake shore/telescopic pole / 0.5–1.5 m below surface	Summer 2009
Fault lakes	Opasatica		Summer 2013
	Rémigny
	Témiscamingue	Boat/telescopic pole/0.5–1.5 m below surface
	Tee	Boat/plastic tubing equipped with an inertial valve/0.5–1.5 m below surface	Autumn 2013
Lakes on clay	Macamic	Sampling at the outlet of the lake/telescopic pole/0.5–1.5 m below surface	Summer 2013
	Abitibi
Rivers	Outaouais	Boat/telescopic pole/central part of the main channel/0.5–1.5 m below surface
	Duparquet	Sampling near the central part of the main channel/telescopic pole or bailer/0.5–1.5 m below surface
	La Sarre
	From Mance lake
	Between Beauchastel and Montbeillard lakes
	Between Routhier and Vallet lakes
	Between Rondelet and Des Quinze lakes
	Kinojevis
Small streams	–	Sampling near the central part of the main channel/telescopic pole or bailer/0.5–1.5 m below surface	Summer/autumn 2013

All surface water samples were filtered immediately after collection using 0.45-μm Waterra FHT. The samples collected for stable isotope analyses (δ²H-δ¹⁸O) were stored in 30–60-mL HDPE bottles and kept at 4°C until analysis.

Springs

Groundwater springs were sampled between July 2009 and September 2011 through the course of a previous study (Castelli 2012). These samples are associated with the Vaudray-Joannès, Launay and Saint-Mathieu-Berry eskers and with the Harricana Moraine. Other spring samples were collected during the summer of 2013. These correspond to springs that are equipped for drinking water supply. All spring samples were filtered immediately after collection using 0.45-μm Waterra FHT carthridges. The in situ parameters (pH, dissolved oxygen, specific electrical conductivity [SEC], redox potential and temperature) were measured on site using a multi-parameter probe. Further details about the sampling methodology can be found in Castelli (2012). The samples collected for stable isotope analyses (δ²H-δ¹⁸O) were stored in 30–60-mL HDPE bottles. The samples collected for tritium (³H) analyses were stored in 1-L HDPE bottles. Sampling bottles for stable isotope analyses (δ²H-δ¹⁸O) and tritium (³H) analyses were filled to the top with no head space and kept at 4°C until analysis.

Groundwater

Municipal, private and observation wells were sampled during this study. The targeted wells were commonly identified using the Quebec Ministry of the Environment water well databases, including the Hydrogeological Information System (HIS). Complementary information about the characteristics of the sampled wells was obtained from the private owners or from municipal employees. Well logs were used to document the characteristics of the sampled observation wells. Two conditions had to be met prior to groundwater sample collection: the infrastructures had to be purged of a minimal volume of water and the monitored in situ parameters had to reach predefined stability criteria. Water was therefore systematically purged from domestic and observation wells before sample collection. This procedure allows draining the stagnant water from the infrastructure to obtain a water sample that represents groundwater from the aquifer. For private wells, it was assumed that water was withdrawn on a regular basis and that the minimal purge could be limited to the volume of the pressurized reservoirs. The purging was done from taps supplying untreated water. The minimum purging volume for observation wells corresponded to 3 times the total volume of water stored within the well and its gravel pack. In municipal wells, it was assumed that water is withdrawn at a sufficient rate to prevent stagnation within infrastructures and that no minimal purge was required. During the purge, in situ parameters (pH, dissolved oxygen, SEC, redox potential and temperature) were measured using a multi-parameter probe. The samples were retrieved only after all of the in situ parameters had reached the stability criteria reported in Table 4 for three consecutive measurements separated by 5-minute intervals. The samples collected for stable isotope analyses (δ²H-δ¹⁸O) were stored in 30–60-mL HDPE bottles and kept at 4°C until analysis. The samples collected for tritium (³H) analyses were stored in 1-L HDPE bottles. Sampling bottles for tritium (³H) analyses were filled under water in a bucket and capped under water with no head space, following a filling procedure described by the United States Geological Survey (USGS), and kept at 4°C until analysis.

Table 4. Stability criteria for groundwater sample collection. Groundwater samples were collected only when the stability criteria were reached for three consecutive measurements separated by 5-minute intervals.

In situ parameters	Stability criteria
Temperature (°C)	± 0.1
pH (pH unit)	± 0.1
Specific electrical conductivity (μS/cm)	± 2%
Dissolved oxygen (mg/L)	± 0.2
Redox potential (mV)	± 10
Turbidity (visual observation)	Constant color/opacity

Analytical procedures

Water samples were analyzed for δ²H and δ¹⁸O at the stable isotopes laboratory of the Geotop-Université du Québec à Montréal (UQAM) research centre in Montreal, Canada. Measurements were made using a dual inlet Micromass Isoprime™ isotope ratio mass spectrometer coupled to an Aquaprep™ system. For oxygen and hydrogen isotopic analyses, 200 μL of water was transferred into septum vials and equilibrated in a heated rack with a known volume of CO₂ (H₂ in the case of hydrogen analyses, with a platinum catalyst). The isotopic compositions of samples were corrected using internal reference waters (δ¹⁸O = −6.71‰, −13.98‰ and −20.31‰; δ²H = −51.0‰, −99.0‰ and −155.4‰) calibrated on the Vienna Standard Mean Ocean Water–Standard Light Antarctic Precipitation (VSMOW-SLAP) scale (Coplen 1996). Reference waters were systematically analyzed between series of samples in order to check for instrumental stability. The analytical uncertainty is ≤ 1‰ for δ²H and ≤ 0.05‰ for δ¹⁸O at the 1σ level. Values are reported in permil units (‰) against the Vienna Standard Mean Ocean Water standard (VSMOW). Water samples for tritium (³H) analyses were measured by liquid scintillation counting (LSC) at the Environmental Isotope Laboratory of the University of Waterloo, Canada. The detection limit for enriched samples is 0.8 tritium units (TU).

Data organization and statistical analyzes

In order to support the interpretation of results, the isotopic data are organized according to the main hydrological components and aquifers that were presented previously (Table 1; Figure 2). Overall, eight distinct groups are identified for segregating the isotopic data. The latter include (with numbers corresponding to those reported in Figure 2) samples corresponding to (1) precipitation, (2) the snowpack, (3) surface waters, (4) springs, (5) unconfined granular aquifers, (6) confined granular aquifers, (7) unconfined fractured rock aquifers and (8) confined fractured rock aquifers. The basic statistical analyses conducted on the data set and related plots were performed using the Statistica version 6 (StatSoft, Inc. 2004) and Microsoft Excel software.

Results and discussion

The entire data set associated with stable isotope analyses is shown in Figure 3a, whereas Figure 3b allows a better representation of surface and groundwater data. Figure 3c presents the average values of groundwater and surface water groups in the same range as Figure 3b. As a complement, Figure 4 presents boxplots showing the non-outlier range, 25th and 75th percentiles and median SEC, δ¹⁸O and ³H values for selected groups of data. The non-outlier range is defined as maximum and minimum values within 1.5 times the interquartile range. It was chosen not to show the data outside of the non-outlier range in the boxplots (Figure 4) in order to allow a better representation of the data. Nevertheless, the entire data set is shown in Figure 3. Results associated with three of the main components of the water cycle (namely precipitation, surface waters and groundwater) are first discussed separately to address specific topics, and jointly interpreted afterward to propose a regional-scale conceptual model of groundwater flow.

Figure 3. Isotopic composition of precipitation, surface waters, springs and groundwater. All data are shown in (a) whereas figure (b) provides a zoom on the surface and groundwater data. Average values are provided in figure (c). The flux-weighted average for precipitation is based on the three stations (Sainte-Hélène-de-Mancebourg, Amos, Bearn) where precipitation samples were collected between 2013 and 2015. The values correspond to the mean of the flux-weighted averages calculated for each station. The arithmetic average calculated for all precipitation samples is also shown. The error bars in (c) correspond to the standard deviation calculated for each sample series. The local meteoric water line (LMWL) and local evaporation line (LEL) derived from the data are also shown. Isotopic values are reported in permil units (‰) against the Vienna Standard Mean Ocean Water standard (VSMOW).

Figure 4. δ¹⁸O (a, b), ³H (c) and SEC (d) boxplots for precipitation, surface waters, springs and groundwater. Outlier data are excluded; the non-outlier range is defined as maximum and minimum values within 1.5 times the interquartile range. The available δ¹⁸O data are shown in (a) whereas a clearer view of surface and groundwater data is provided in (b). The classification reported in (c) is adapted from Clark and Fritz (1997). UG: unconfined granular aquifers; UFR: unconfined fractured rock aquifers; CG: confined granular aquifers; CFR: confined fractured rock aquifers; SEC: specific electrical conductivity. In (a) and (d), the samples collected from the snowpack and at the precipitation monitoring stations are shown separately. δ18O values (a,b) are reported in permil units (‰) against the Vienna Standard Mean Ocean Water standard (VSMOW). δ¹⁸O values (a,b) are reported in permil units (‰) against VSMOW.

The isotopic composition of precipitation: air masses and seasonal patterns

The δ²H and δ¹⁸O values of precipitation ranged from −217.9‰ to −29.9‰ and −28.4‰ to −4.8‰, respectively (Figure 3a). These data are used to define the local meteoric water line, calculated using precipitation (n = 164) and snowpack data (n = 31) (LMWL; δ²H = 7.84 δ¹⁸O + 9.20 ‰ [n = 195, r² = 0.99]; Figure 3). The LMWL is relatively close to Craig’s global meteoric water line (GMWL; δ²H = 8 δ¹⁸O + 10 ‰ Standard Mean Ocean Water (SMOW); Craig 1961). In terms of temporal variations, the monthly precipitation data revealed the expected temperature-dependent isotopic cycle, with heavy isotope depletion during the colder period and enrichment during the warmer period (Figure 5). Such a pattern is typical of boreal regions with summer rains enriched in heavy isotopes and winter precipitation depleted in heavy isotopes in response to strong seasonal variations in temperature (Clark and Fritz 1997). Overall, the temporal variations in the isotopic composition of precipitation reflect the marked seasonality of the study area.

Figure 5. Seasonal variations in precipitation amounts and δ¹⁸O composition. See Figure 1 for the location of monitoring stations. The histogram showing precipitation amounts is based on data from July 2013 to July 2015 as recorded at the monitoring stations. Note: SHM: Sainte-Hélène-de-Mancebourg. δ¹⁸O values are reported in permil units (‰) against the Vienna Standard Mean Ocean Water standard (VSMOW).

The weighted mean annual isotopic compositions of precipitation at stations 1, 3 and 4 are shown in Table 5. The data do not show a clear relationship between precipitation δ²H-δ¹⁸O and latitude. The absence of a latitudinal gradient in precipitation isotopic composition could reflect the superimposed influence of the three major air masses that supply moisture over the study area, namely (1) the Arctic air mass (from the north), (2) the tropical vapor (from the south) and (3) the westerlies (Bryson 1966; Fritz et al. 1987). Given the relatively flat landscape of the study area, the effect of altitude on precipitation δ²H-δ¹⁸O values is probably limited. Finally, evaporation processes might also influence the isotopic composition of precipitation at the local scale, near large water bodies such as the Témiscamingue, Simard and Abitibi lakes (Figure 1). Further data would be required to better explain spatial variations in the isotopic composition of precipitation at the scale of the study area.

Table 5. Amount weighted average isotopic composition of precipitation at the four sampling stations (2013–2015).

Stations	n	δ^2H (‰ VSMOW)	δ^18O (‰ VSMOW)
Station 1 (Béarn)	24	−91.9	−13.0
Station 2 (Montbeillard)	13	NA	NA
Station 3 (Sainte-Hélène-de-Mancebourg)	24	−100.9	−14.2
Station 4 (Amos)	24	−84.5	−12.0

Notes: Refer to Figure 1 for station locations; n: number of samples.

NA: averages values are not provided owing to a lack of data from November to April.

VSMOW: Vienna Standard Mean Ocean Water standard.

Samples collected from the snowpack in March 2014 (prior to the onset of spring snowmelt) revealed δ²H and δ¹⁸O values ranging from −184.4‰ to −151.6‰ and –24.5‰ to –20.6‰, respectively (Figure 3). Ice layers have been found in some places in the snowpack during sample collection. It is recognized that isotopic fractionation occurs during snowmelt (Taylor et al. 2001; Laudon et al. 2002). This suggests that the different layers identified within the sampled snowpack (which reflect melting/refreezing events) are most likely characterized by variable isotopic compositions. Nevertheless, here, the bulk samples (i.e. composite samples of the entire thickness of snow) collected from the snowpack prior to the onset of the snowmelt event plot uniformly on the left side of the LMWL (Figure 3), with values similar to those measured at the precipitation monitoring stations. This suggests that the bulk of the snowpack preserved an isotopic composition that is similar to that of precipitation prior to the onset of the spring snowmelt event. Overall, the isotopic compositions of samples retrieved from the snowpack do not define a clear latitudinal gradient at the regional scale.

The isotopic composition of surface waters: inflows, evaporative enrichment and mixing

The regional data set reveals a significant variability in surface water δ²H and δ¹⁸O, with values ranging from −100.8‰ to −68.9‰ and −13.7‰ to −7.6‰, respectively (for regional average values of −80.9‰ [δ²H] and −10.7‰ [δ¹⁸O]). The isotopic composition of surface waters can present both spatial and temporal variations owing to the superimposed influence of precipitation, groundwater inflow, evaporation and tributary mixing processes (Yi et al. 2010). Temporal variations are commonly observed in regions characterized by a marked seasonality where heavy isotope depletion following snowmelt and enrichment caused by evaporation during the ice-off period have been reported (e.g. see Telmer and Veizer 2000; Yi et al. 2010). Here, most of the samples were collected in summer (Table 2), and the data plot below the LMWL, consistent with the hypothesis that surface waters have undergone heavy isotope enrichment owing to evaporation. A surface local evaporation line (LEL): δ²H = 3.80 δ¹⁸O – 40.22 ‰ (n = 99, r² = 0.65; Figure 3a) is evaluated from the data. The evaporation over inflow ratio (E/I), as described by Gibson and Edwards (2002), corresponds to the water balance of a surface water body and relates to the combined effects of evaporative heavy isotope enrichment and dilution by inflowing waters. The E/I ratio is used here for evaluating the hydrological processes responsible for the distribution of surface water data plotting below the LMWL in δ²H vs. δ¹⁸O graphs (Figure 3). The Hydrocalculator application from Skrzypek et al. (2015) was used to calculate E/I ratios under steady-state conditions (shown over the LEL in Figure 6a). This application solves a modified form of the Craig and Gordon (1965) model:(2)

Figure 6. Interpretation of the isotopic compositions of surface waters. Data in (a) represent the local evaporation line evaluated using the 99 available surface water samples. The Evaporation over Inflow (E/I) ratios are shown for points plotting on the Local Evaporation Line (LEL) at 0.5‰ intervals along the x-axis. These E/I ratios are shown for three scenarios (see text for details). Further details with respect to the types of surface water samples are provided in (b) and in Table 3. Isotopic values are reported in permil units (‰) against the Vienna Standard Mean Ocean Water standard (VSMOW).

where δ_L corresponds to the isotopic composition of the sampled surface water, is the initial isotopic composition of water and δ^* is the limiting isotopic enrichment factor (a function of the isotopic fractionation factor and atmospheric isotopic composition and relative humidity). Similar approaches are used to evaluate the E/I ratios of rivers (e.g. Telmer and Veizer 2000) and lakes (e.g. Gibson and Edwards 2002). Here, theoretical calculations are proposed for points falling directly on the LEL evaluated from the data set, the latter corresponding to the δ_L term of Equation (2)(2) . The initial isotopic composition of water (δ_P) is assumed to correspond to the intersection between the LMWL and the LEL. The isotopic composition of atmospheric moisture (δ_A), which is needed to evaluate in Equation (2)(2) , is evaluated according to three scenarios:

(1)	δA is in equilibrium with the flux-weighted average isotopic composition of precipitation between the months of May and August (δA = −141.21 ‰; −19.45 ‰);
(2)	δA is in equilibrium with the arithmetic average isotopic composition of precipitation between the months of May and August (δA = −137.42 ‰; −18.91 ‰);
(3)	δA is in equilibrium with δP (δA = −156.03 ‰; −21.63 ‰).

A temperature of 22°C and a relative humidity of 70% are assumed, consistent with the conditions encountered during the warm ice-off period within the study region. The isotope fractionation factors are calculated from the Horita and Wesolowski (1994) equations.

Depending on the model parameters used, E/I ratios ranging between 0 and 36% are evaluated from the regional data set (Figure 6a). Based on this theoretical framework, surface water samples plotting closer to the LMWL most likely present the lowest E/I ratios. These surface waters might be affected by significant evaporation losses, but inflows (mostly precipitation and groundwater, both plotting along the LMWL) are sufficient to buffer the heavy isotope enrichment resulting from evaporation. On the contrary, the samples that are most enriched in heavy isotopes (i.e. plotting farther along the right side of the LEL) most likely present the highest E/I ratios. Within the regional hydrogeological framework, it is most likely that evaporation mainly occurs within lakes and stagnant water bodies, whereas dynamic in-stream processes mainly reflect surface water/groundwater interactions and tributary mixing.

One distinctive feature of the regional data set is that most of the kettle lakes plot on the right side of the LEL, suggesting overall high E/I ratios in comparison to other surface waters (Figure 6b). Recent studies by Isokangas et al. (2015) and Arnoux et al. (2017) provided novel methods to estimate the groundwater dependency of such lakes based on E/I ratios calculated from stable isotope data, among other things. Here, quantifying the groundwater dependency of lakes is precluded by uncertainties related to volumes, surface water inflows and potential heterogeneities caused by thermal stratification and temporal variations. Nevertheless, the distinctive isotopic composition of kettle lakes stresses the need to conduct further studies aimed at quantifying their water balance to better anticipate their sensitivity to climate change and human pressure.

The isotopic composition of groundwater and springs: recharge, discharge, mixing and regional flowpaths

The stable isotope analyses (δ²H-δ¹⁸O) conducted on 351 groundwater and spring samples are shown in Figures 3, 4a and 4b. The data associated with groundwater samples are represented according to the four main types of regional aquifers, as shown in Figure 2 (unconfined granular aquifers; unconfined fractured rock aquifers; confined granular aquifers; confined fractured rock aquifers). The δ²H and δ¹⁸O values of groundwater samples ranged from −107.8‰ to −71.9‰ and −15.2‰ to −9.4‰, respectively, for a regional average of −13.3‰ (δ¹⁸O) and −93.4‰ (δ²H). The δ²H and δ¹⁸O values of spring samples ranged from −106.3‰ to −75.1 ‰ and −14.4‰ to −9.7‰, respectively, for a regional average of −13.8‰ (δ¹⁸O) and –98.1‰ (δ²H). Previous work by Castelli (2012) revealed that the isotopic composition of springs is fairly stable (1) throughout the ice-off season and (2) from year to year. Similarly, it is assumed here that the isotopic composition of groundwater is relatively constant in time. The focus is therefore set on spatial rather than temporal variations in isotopic compositions. The data associated with groundwater and springs are evenly distributed around the LMWL. This suggests that (1) infiltration and recharge processes are not significantly affected by evaporation losses and (2) aquifers are not significantly affected by inflows from evaporated surface waters. This is consistent with the regional hydrogeological framework because (1) the main recharge areas consist of coarse-grained glaciofluvial formations where high infiltration rates prevent surface pounding (and associated evaporation) and (2) physical measurements suggest that groundwater generally discharges into surface waters at the regional scale (Cloutier et al. 2013, 2015). Although water losses by plant transpiration prior to groundwater recharge might be significant, the available data set does not allow the quantification of this process. Results also support the hypothesis that regional surface waters – groundwater interactions are essentially restricted to groundwater discharging into surface waters, and not the inverse (otherwise, some groundwater samples would most likely be plotting significantly below the LMWL owing to mixing with evaporated surface waters). The data do not show systematic latitudinal trends in the stable isotope composition of groundwater at the regional scale. This is most likely caused by the complex mixing of recharge originating from the main atmospheric air masses supplying precipitation over the study region.

The data associated with the four main types of aquifers present a significant dispersion and overlap within Figures 3 and 4b. Nevertheless, two main observations arise from the distribution of data within the δ²H-δ¹⁸O plots:

(1)	The data associated with unconfined granular aquifers present a greater dispersion in comparison to the other types of aquifers (Figure 4b);
(2)	Confined granular aquifers and fractured rock aquifers (both confined and unconfined) tend to plot farther on the right side of LMWL in comparison to springs and unconfined granular aquifers (Figure 3c).

These observations are interpreted in an effort to further constrain the understanding of regional flowpaths within the hydrogeological environment conceptually represented in Figure 2. The unconfined granular aquifers represent the main regional recharge zones where rain and melted snow preferentially infiltrate at different periods of the year, depending on local climatic and hydrological conditions (Cloutier et al. 2016). The proposed interpretation of groundwater stable isotope data is that unconfined granular aquifers present a greater δ²H-δ¹⁸O variability owing to a shorter residence time of water within shallow and dynamic systems such as eskers and moraines. The results associated with ³H analyses (Figure 4c) further support this interpretation. Overall, the tritium data suggest that water samples from springs, unconfined granular aquifers and unconfined fractured rock aquifers correspond to modern waters (with values ranging between approximately 9 and 17 TU, consistent with the subdivisions proposed by Clark and Fritz 1997). Such results suggest a relatively short residence time of water in these shallow components. In contrast, samples collected within confined aquifers (both granular and fractured) present tritium contents below 0.8 TU (consistent with submodern values according to the subdivisions proposed by Clark and Fritz 1997). Such results suggest that within the region, confined aquifers generally contain older groundwater in comparison to unconfined aquifers.

The sampled springs are all located on the flanks of unconfined granular aquifers associated with glaciofluvial formations (eskers/moraines). Their heavy isotope depletion (Figure 3c) with respect to most groundwater samples suggests that superficial flow systems associated with springs are submitted to a greater influence of lighter precipitation. As previously discussed, within the study area, precipitations depict a clear seasonal pattern, with lighter values recorded during the cold period (Figure 5). Under such conditions, one likely explanation is that a significant proportion of the groundwater recharge associated with snowmelt contributes to spring discharge. This stands as a plausible explanation because the springs located on the flanks of eskers and moraines essentially consist in ‘overflows’ allowing groundwater to discharge from unconfined granular aquifers towards the surface water network. From a physical perspective, it seems likely that ‘overflow’ conditions would be favored following the thawing period, as groundwater levels are generally rising within the eskers and moraines of the region during this period (Cloutier et al. 2013, 2015). During the drier summer period, recharge rates are likely reduced, therefore limiting the supply of heavier water associated with summer precipitation to shallow flow systems. In fall, at the end of the warmer period, when the uptake of water by vegetation is reduced and precipitation is increased, recharge rates most likely increase, further supplying water that is depleted in heavy isotopes to the shallow flow systems associated with springs. The loss of lighter isotopes resulting from groundwater discharge within springs could also partly explain why the confined aquifers, located farther downstream within the flow systems, are slightly enriched in heavy isotopes with respect to springs and shallow unconfined aquifers (Figures 3c and 4b). One key regional observation discussed by Cloutier et al. (2016) is that extensive peatlands are developed on the flanks of several eskers and moraines of the region. The contacts between eskers/moraines and these extensive peatlands are interpreted as diffuse groundwater seepage zones (Nadeau et al. 2015). The process related to the loss of lighter isotopes within punctual springs might therefore be generalized to most of the unconfined aquifers associated with extensive peatlands and diffuse groundwater seepage zones. In addition, it seems most likely that the areas associated with unconfined fractured rock aquifers are less prone to groundwater recharge during snowmelt in comparison to the areas associated with glaciofluvial formations (i.e. unconfined granular aquifers). This hypothesis is based on the contrast in hydraulic conductivity (K) between the fractured bedrock (lower K) and the eskers and moraines (higher K; e.g. see Cloutier et al. 2015). The differences in the isotopic composition of unconfined granular aquifers and springs compared to that of confined granular aquifers and fractured rock aquifers (both confined and unconfined) could therefore reflect the relative contribution of snowmelt-induced recharge for the different types of aquifers.

The interpretations discussed above are further illustrated in Figure 7a, where the SEC of springs and groundwater is plotted against the corresponding δ¹⁸O values. The data presented therein reveal that confined fractured rock aquifers generally present the highest SEC values, whereas the unconfined granular aquifers present the lowest SEC values. This is also depicted in Figure 4d where a gradual increase in SEC is observed from unconfined towards confined aquifers. SEC can thus be used as a proxy for groundwater residence time and geochemical evolution, higher SEC generally implying more evolved groundwater. This pattern most likely reflects the combined effects of (1) increased geochemical interactions between groundwater and the geological environment along regional flowpaths and (2) recharge of precipitation having very low dissolved solid concentration (thus quite low SEC) in unconfined granular aquifer areas. The mixing model shown in Figure 7a further supports the interpretations proposed above. Three end-members are proposed to explain the variability in both δ¹⁸O and SEC values:

(1)	End-member ‘A’ presents low SEC and δ18O; it is identified as representative of recent snowmelt-induced recharge. The lowest non-outlier δ18O value recorded in unconfined granular aquifers (lower end of the boxplot shown in Figure 4b) is used to characterize this end-member, whereas the SEC value is set to 0.05 mS/cm, corresponding to the lowest non-outlier SEC value recorded in unconfined granular aquifers (lower end of the boxplot shown in Figure 4d);
(2)	End-member ‘B’ presents low SEC and high δ18O close to −12‰; it is identified as representative of fall recharge affected by relatively high temperature rain. The highest non-outlier δ18O value recorded in unconfined granular aquifers (higher end of the boxplot shown in Figure 4b) is used to characterize this end-member, whereas the SEC value is set to 0.05 mS/cm, as previously defined for end-member A;
(3)	End-member ‘C’ presents high SEC and δ18O; it is identified as representative of well-mixed water (from various recharge zones or various recharge periods) that is geochemically evolved owing to prolonged water–rock interactions along flowpaths, thus leading to a higher SEC. The median δ18O isotopic composition of confined fractured rock aquifers is used to constrain this end-member, and the SEC is set to 1.8 mS/cm (the highest recorded value in the data set; Figure 4d).

The dashed lines in Figure 7a and b represent mixing lines between the three end-members, at 25% intervals. In this mixing model, the SEC is considered linearly related to the geochemical level of evolution of groundwater, which is most probably not linearly related to the groundwater residence time.

Figure 7. Specific electrical conductivity (SEC) vs δ¹⁸O of springs, surface waters and groundwater. Three end-members are shown: (A) corresponds to groundwater recharged during snowmelt, (B) corresponds to groundwater that is less impacted by snowmelt-induced recharge and (C) corresponds to geochemically evolved groundwater found farther along regional flowpaths. See text for further details. The six main outliers (as described in the text) are not shown in this figure. Figure (a) shows the distribution of data whereas the corresponding averages are shown in (b). δ¹⁸O values are reported in permil units (‰) against the Vienna Standard Mean Ocean Water standard (VSMOW).

Six main outliers are identified among the groundwater/springs data set in Figures 3b. These six points plot below the LMWL in Figure 3, suggesting an evaporative enrichment that is not observed in the other groundwater samples. Five of these points are associated with samples collected in the immediate vicinity of lakes, whereas the sixth one was collected within a mining site. It is proposed that the five samples collected near lakes and presenting heavy isotope enrichment and low SEC values are directly affected by the supply of evaporated and low-SEC lake waters. The wells characterized by low SEC and low δ¹⁸O values are likely affected by nearby snowmelt-induced groundwater recharge, whereas those with low SEC and significantly higher δ¹⁸O values are likely affected by mixing with surface waters. In that sense, within the study region, stable isotope data could be used as a screening tool to detect groundwater wells that are affected by surface waters. Such wells would deserve specific attention in groundwater monitoring programmes because they might be more vulnerable to microbial contamination originating from surface waters. Finally, the sixth outlier point was collected in a groundwater well located within a mining site; its isotopic composition and SEC value most likely reflect anthropogenic effects.

Local and regional systems of groundwater flow: insights from water stable isotopes

Tóth (1963) proposed a conceptual model of groundwater flow in two-dimensional cross sections involving local, intermediate and regional groundwater flow sub-systems. Such conceptual models constitute a central unifying concept to describe groundwater recharge, flow and discharge processes at various scales and provide a framework that is of upmost importance to interpret regional hydrogeochemical data sets. Tóth (1999) further developed his conceptual model by presenting groundwater as a geologic agent interacting with the environment and inheriting a geochemical composition that preserves an archive of hydrogeochemical processes and systematically evolves along flowpaths. Following Tóth’s pioneer work, the different flow systems of the study area are interpreted with a focus on the isotopes of the water molecule itself, rather than other dissolved species. The working hypothesis is that within the considered spatiotemporal domain, the isotopic composition (δ²H-δ¹⁸O) of water is not significantly affected by water–rock interactions but rather mainly reflects the superimposed effects of precipitation isotopic composition, evaporation, recharge, mixing and discharge processes that vary in time and space. The average δ²H-δ¹⁸O compositions of precipitation (including snow samples), surface waters, springs and groundwater are shown in Figure 8, which consists in a more detailed representation of the regional hydrogeological conditions that were initially conceptualized in Figure 2.

Figure 8. Isotopic variations within local, intermediate and regional flow systems. The range of measured isotopic compositions is shown for the sampled hydrological components and hydrogeological units, as identified in Tables 1 and 2. Areas of preferential snowmelt-induced recharge are identified. Flowpaths at various scales are shown. (R) Recharge area; (E) evaporation; (T) transpiration; (GW/SW) groundwater/surface water interactions (herein mainly groundwater exfiltration). The numbers (1 to 8) correspond to the sampled hydrological components and hydrogeological units, as identified in Table 1.

The conceptual model illustrates precipitation showing the greatest isotopic variability and recharge occurring mainly within shallow unconfined aquifers, where the variability in the isotopic composition of water is attenuated with respect to precipitation due to mixing. Part of the water flowing within these unconfined aquifers discharges towards surface waters through springs and diffuse seepage zones that are generally depleted in heavy isotopes with respect to the average composition of groundwater. Surface waters are subsequently submitted to evaporation (causing heavy isotope enrichments) and to further dilution by precipitation and groundwater discharge. The remaining groundwater flows from unconfined recharge areas towards confined aquifers, allowing two key geochemical processes to occur along flowpaths:

(1)

Water–rock interaction: groundwater becomes increasingly concentrated in dissolved ions owing to prolonged water–rock interactions, as shown by the gradual increase in SEC among the various hydrological components, from precipitation towards confined fractured rock aquifers (Figure 4d). This interpretation is further supported by the tritium analyses (Figure 4c), which suggest that confined aquifers are characterized by longer water residence times compared to springs and unconfined aquifers;

(2)

Mixing: The isotopic composition of groundwater evolves towards a regional average owing to the mixing of water originating from various recharge zones or various recharge periods. This is illustrated in Figure 4b and Figure 7a, which suggest that mixing processes tend to damp the variability recorded in δ18O values within confined aquifers. In other terms, the greater the proportion of unconfined aquifers in the vicinity on a given site, the greater the potential variability in the stable isotope composition of the groundwater.

Clark and Fritz (1997) proposed a schematic representation of the attenuation of seasonal isotope variations in recharge waters within the unsaturated zone. These authors reported that temporal isotopic variations around an average value are gradually damped with increasing depth as a function of flowpath lengths and residence times. A comparable interpretation is proposed here to explain the evolution of groundwater SEC and δ²H-δ¹⁸O along local, intermediate and regional flow systems (Figures 7 and 8). On one hand, since recharge processes are primarily associated with sectors of unconfined aquifers and partly controlled by hydrogeological conditions that are highly heterogeneous at the regional scale, a significant variability is recorded in unconfined aquifers’ isotopic compositions. Overall, the unconfined granular aquifers associated with coarse-grained glaciofluvial formations (eskers and moraines) are more prone to recharge during snowmelt and therefore present low δ²H-δ¹⁸O values in comparison to the other types of aquifers. Similarly, the springs found on the flanks of eskers and moraines are depleted in heavy isotopes. On the other hand, recharge areas associated with unconfined bedrock aquifers are more prone to runoff generation during the intense snowmelt event, and are therefore generally enriched in heavy isotopes in comparison to unconfined granular aquifers. Finally, the confined aquifers located farther along groundwater flowpaths are fed by numerous interconnected unconfined aquifer areas that are heterogeneously distributed over the region. The result of the mixing processes occurring along flowpaths is a gradual evolution towards a regional average value within confined aquifers. Based on this interpretation, groundwater samples from confined aquifers that present an isotopic composition that significantly differs from the regional average most likely represent sites where the preferential influence of a nearby recharge area are recorded.

Conclusions

The isotopic composition of water (δ²H-δ¹⁸O and ³H) in precipitation, surface waters, groundwater and springs was used to improve the understanding of regional flow systems in the southern portion of the Barlow-Ojibway Clay Belt. The data allowed the identification of the key hydrological processes driving spatiotemporal isotopic variations and aided in deciphering the interaction pathways between the various components of the water cycle, from atmospheric moisture to groundwater discharge:

(1)	Regional precipitation depicts a clear seasonal pattern with heavy isotope depletion recorded in winter and heavy isotope enrichment found in summer. The data did not reveal clear spatial trends in precipitation δ2H-δ18O, most likely in response to the complex mixing of inputs from the main air masses supplying moisture to the study region. Overall, precipitation data allowed the definition of the local meteoric water line (LMWL: δ2H = 7.84 δ18O + 9.20 ‰; n = 164; r2 = 0.99);
(2)	Surface waters plot below the LMWL in a δ2H-δ18O graph owing to evaporation losses. Calculations conducted using a modified form of the Craig and Gordon (1965) model with the Hydrocalculator application (Skrzypek et al. 2015) provided an estimation of evaporation over inflow (E/I) ratios ranging between 0 and 36% at the regional scale;
(3)	Groundwater samples plot on the LMWL and are evenly distributed around the average isotopic composition of precipitation in a δ2H-δ18O graph, suggesting that evaporation processes do not significantly affect regional recharge waters. The δ2H-δ18O of groundwater in shallow unconfined aquifers presents a greater variability in comparison to that of confined aquifers owing to mixing processes that tend to buffer the isotopic variability along the regional groundwater flowpaths. This buffering of the isotopic composition variability is accompanied by increasing values of the SEC of groundwater from unconfined to confined aquifers, most likely in response to increasing water–rock interaction timescales, consistent with generally lower tritium contents in confined aquifers;
(4)	Unconfined granular aquifers and associated springs are generally slightly depleted in heavy isotopes in comparison to other groundwater samples, suggesting that these shallow flow systems are preferentially recharged during spring snowmelt within the region.

Overall, the observations and interpretations of this study support an improvement of a conceptual hydrogeological model of local and regional hydrological systems and flowpaths that is consistent with the geological framework. This conceptual model illustrates how water isotopes can be used to constrain physically based interpretations of local, intermediate and regional flow systems such as those inspired by the pioneering work of Tóth (1963), even in a highly heterogeneous region. The buffering of isotopic variability in groundwater from groundwater recharge areas located in unconfined aquifers towards groundwater discharge areas and confined aquifers located farther along regional flowpaths provides a framework for deciphering the processes associated with the geochemical evolution of groundwater, as a complement to understandings based on dissolved major ions and trace elements.

Acknowledgements

The authors would like to acknowledge the financial contribution of the Quebec Ministry of the Environment (Ministère du Développement durable, de l’Environnement et de la Lutte contre les changements climatiques) through the Programme d’Acquisition des Connaissances sur les Eaux souterraines (PACES). We also acknowledge the financial contribution of regional partners involved in the PACES programme, including the Regional County Municipalities (Abitibi, Vallée-de-l’Or, Abitibi-Ouest, Ville de Rouyn-Noranda, Témiscamingue) and the Regional Conference of Elected Officials of Abitibi-Témiscamingue. The authors acknowledge the Fondation of the University of Quebec in Abitibi-Témiscamingue (FUQAT) for scholarships to the project of Nathalie Rey. We thank Magalie Roy for geomatics support. The work of students involved in the project as field assistants, and the collaboration of the population of Abitibi-Témiscamingue giving site access, are greatly appreciated. The authors kindly thank two anonymous reviewers for constructive comments that greatly contributed to improve the original manuscript.