Benefits and limitations of using isotope-derived groundwater travel times and major ion chemistry to validate a regional groundwater flow model: example from the Centre-du-Québec region, Canada

Abstract

Understanding groundwater dynamics at the regional scale (> 100 km) is essential to the development of sustainable water management regulations. Groundwater flow models are increasingly used to support these strategies. However, in order to be reliable, these models need to be calibrated and validated. The objective of this work is to evaluate the benefits and the limitations of using isotope-derived groundwater travel times and major ion chemistry to validate a regional-scale groundwater flow model in the humid continental climate of southern Québec (Canada). A three-dimensional regional-scale steady-state groundwater model was created using MODFLOW for the fractured bedrock aquifer of the Centre-du-Québec region (Québec, Canada), using data acquired during recent aquifer characterization projects. The model covers an area of 7452 km², from the unconfined Appalachian Mountains to the confined St. Lawrence Platform. Groundwater travel times were simulated for 211 wells using particle tracking. The groundwater flow model was calibrated using 11,775 regionally distributed heads and 15 baseflow values. The model was validated using 23 ³H/³He residence time (3 to 60 years), 17 ¹⁴C residence time (226 to 22,600 years), and the major ion compositions from 211 wells. Results indicate that the model is able to satisfactorily simulate ³H/³He isotopic residence time, while ¹⁴C isotopic residence times are generally underestimated. These results suggest substantial mixing between groundwater recharged during the last deglaciation and recently recharged water. Regional groundwater flow is limited or absent, and most of the recharge discharges to the river network as baseflow. The analysis of travel times indicates a statistically distinct mean travel time for the different groundwater types. Median travel time is 68 years for recently recharged groundwater (Ca-HCO₃), 274 years for semi-confined groundwater (Na-HCO₃), and 738 years for confined groundwater (Na-Cl). This confirms that groundwater chemistry is a broad indicator of groundwater travel time.

La compréhension de la dynamique régionale (> 100 km) de l’eau souterraine est essentielle au développement d’une règlementation orientée vers le développement durable de cette ressource. Les modèles d’écoulement de l’eau souterraine sont de plus en plus utilisés pour supporter ces stratégies. Par contre, pour être utilisés à des fins de règlementation, ces modèles doivent être calés et validés. L’objectif de ce travail est d’évaluer les avantages et les limites de l’utilisation de l’âge isotopique de l’eau souterraine et de la géochimie des ions majeurs pour valider un modèle régional de l’écoulement de l’eau souterraine dans le climat continental humide du sud du Québec (Canada). En utilisant les données acquises dans le cadre de projets de caractérisation hydrogéologiques récents, un modèle à trois dimensions régional en régime permanent a été construit avec MODLOW afin de représenter l’aquifère fracturé de la région du Centre-du-Québec (Québec, Canada). Le modèle couvre une superficie de 7 452 km² à partir des zones en nappe libre des montagnes Appalachiennes jusqu’aux secteurs de nappe captive de la plate-forme du Saint-Laurent. Les temps de parcours de l’eau souterraine ont été simulés à 211 puits en utilisant le traçage de particules (MODPATH). Le modèle d’écoulement souterrain a été calé avec 11 775 niveaux distribués dans toute la zone d’étude et avec 15 mesures de débit de bases obtenues à partir de séparation d’hydrogramme ou de mesures manuelles en période d’étiage. Le modèle a ensuite été validé en utilisant les âges isotopiques dérivés de 23 ratios ³H/³He (entre 3 et 60 ans) et de 17 valeurs de ¹⁴C (entre 226 et 22 600 ans), ainsi qu’avec 211 résultats de géochimie des ions majeurs. Les résultats montrent que le modèle est en mesure de simuler de manière satisfaisante les âges ³H/³He tandis que les âges ¹⁴C sont généralement sous-estimés. Ces résultats suggèrent un mélange important entre l’eau rechargée durant la dernière déglaciation et l’eau récemment rechargée. Le flux d’eau souterraine régional est faible ou absent et la grande majorité de la recharge retourne au cours d’eau sous forme de débit de base. L’analyse des temps de parcours moyens indique une différence statistiquement significative (α = 0.01) entre les types d’eau souterraine. Le temps de parcours médian est de 165 ans pour les zones de recharge (Ca-HCO₃), de 746 ans pour les zones semi-captives (Na-HCO₃) et de 7 841 ans pour les zones captives (Na–Cl). Ces résultats montrent bien l’utilité des types d’eau souterraine comme indicateurs généraux des temps de parcours de l’eau souterraine.

Introduction

Understanding groundwater system dynamics at a regional scale (> 100 km) is essential to the development of sustainable water management strategies. Indeed, groundwater systems and aquifer vulnerability and sustainability are increasingly integrated into regional land development plans and climate change impact scenarios (Morris et al. 2003; Reilly et al. 2008; Gleeson et al. 2012). Approaches used to assess the sustainability and vulnerability of groundwater resources include an ecosystemic approach to include all components of the water cycle (Collin & Melloul 2003), the integration of longer time scales (i.e. inter-generational) for aquifer management (Gleeson et al. 2010) and the development of geographic information system (GIS) tools to integrate multi-source data into management tools (Pandey et al. 2011; Przemyslaw et al. 2016). Groundwater flow models are also important tools for hydrogeologists, and have become a regular part of aquifer studies over the past few decades. However, many challenges arise with regards to the development of reliable aquifer models when data from different sources are used (e.g. geology, hydrology, water chemistry, and isotopes), which represent spatial scales ranging from metres to thousands of kilometres, and temporal scales ranging from days to thousands of years (e.g. Scanlon et al. 2002). Increasing the density of aquifer data available at different spatial and temporal scales can reduce these challenges, but even the most extensive aquifer characterization effort will leave numerous issues unresolved due to the inherently heterogeneous nature of aquifers (de Marsily et al. 2005).

Creating regional-scale flow models always implies a certain level of simplification, because the available hydrogeological information is never as dense as the grid resolution of the model. This simplification can create problems of non-uniqueness during model calibration, leading to many different solutions that can reproduce the observed conditions equally well. Calibration based on potentiometric heads alone is considered insufficient to validate the choice of a conceptual flow dynamic model (Voss 2011). Additional field-measured data, such as baseflows, groundwater flow rates, natural and artificial tracer concentrations, or isotopic groundwater residence time, provide a means to better define flow directions. The use of multiple calibration targets is increasingly common in groundwater flow models developed for research purposes (Portniaguine & Solomon 1998; Sanford et al. 2004; Troldborg et al. 2008; Sanford 2011; Turnadge & Smerdon 2014), but is not yet common practice in the industry.

Baseflow data are usually obtained by hydrograph separation using algorithms that make use of the total flow rate time series (e.g. Chapman 1991; Arnold & Allen 1995; Eckhardt 2005). It is generally recognized that these methods tend to overestimate baseflows, especially during high-flow periods, such as spring snowmelt or storm events (Croteau et al. 2010; Rivard et al. 2014). Manual stream flow measurements during low-flow periods are also an effective way to estimate river baseflows, but values obtained in this way represent a single baseflow measurement, which would need to be repeated at different times to provide a reliable estimate of groundwater discharge to streams and rivers. During the model calibration process, estimated baseflows can be compared to river and stream outflows simulated by the groundwater flow model.

Groundwater travel times, or residence time, can be calculated using isotopic dating methods based on tritium (³H) and its product, helium-3 (³H/³He), or using ¹⁴C (e.g. Phillips & Castro 2003). Tritium is the radioactive isotope of hydrogen that decays with a half-life of 12.32 ± 0.02 years to its stable daughter, ³He. Naturally occurring at low levels by cosmic-ray interactions with nitrogen in the atmosphere, tritium was released in large amounts in the 1960s during nuclear weapon testing and introduced into aquifers via recharge. ³H/³He makes it possible to calculate groundwater residence time younger than 70 years (e.g. Tolstikhin & Kamenskiy 1969). The radioactive isotope of carbon, radiocarbon or ¹⁴C, is produced in the upper atmosphere by reaction with nitrogen. Oxidized to ¹⁴CO₂, it enters the hydrological cycle as soil CO₂. Because it decays to stable ¹⁴N in 5730 ± 40 years, groundwater from a few hundred to 35,000 years old can be dated using this method (e.g. Plummer & Glynn 2013).

Because different age tracers date water of different residence time (i.e. younger versus older water), combining tracers is essential to understanding the distribution of groundwater residence time for a given water sample (Suckow 2014) and to identifying mixtures of multiple water masses (e.g. Saby et al. 2016). In theory, groundwater residence time obtained using isotopic tracers can be compared to simulated advective travel times calculated from a numerical groundwater flow model. Although a complete simulation of transport processes, including dispersion and diffusion, might yield travel times that are closer to those estimated using isotopic residence time tracers (e.g. Suckow 2014; Wen et al. 2016), simulated advective travel times can provide an acceptable first estimate. It is generally recognized that groundwater types are representative of flow paths (e.g. Cloutier et al. 2006; Blanchette et al. 2013; Saby et al. 2016) and reflect, to some extent, aquifer vulnerability (Meyzonnat et al. 2016). This comparison is rarely reported in the literature and should be added to model validation when data are available. However, although it is acknowledged that a groundwater flow model might adequately represent the hydrogeological system at a regional scale, the absence of a detailed representation of the local scale heterogeneity might lead to substantial errors in the simulation of groundwater travel times (Larocque et al. 2009). A numerical model unable to significantly differentiate groundwater types or reproduce groundwater travel times reasonably well would indicate a problem with the conceptual model or with the spatial distribution of recharge.

Numerous regional groundwater modelling studies have been carried out in the United States over the last 20 years. The United States Geological Survey (USGS) has funded many projects within the ‘Water Availability’ and the ‘Use Science’ programmes (Reilly et al. 2008), mainly to establish water budgets, identify aquifer stresses, and determine groundwater extraction sustainability. Several projects have also included geochemical and residence time tracers, such as ³H/³He and ¹⁴C, to constrain model calibration (e.g. Sheets et al. 1998; Izbicki et al. 2004; Sanford et al. 2004; Gusyev et al. 2013). However, few regional-scale groundwater flow modelling studies (using water chemistry and residence time tracers) have been reported for the hydrogeological settings of eastern Canada under conditions typical of deglaciation-impacted aquifers and a humid climate. In this area, regional groundwater flow models have been developed to estimate recharge rates or baseflows (Beckers & Frind 2001; Rivard et al. 2014), to assess groundwater sustainability (Meriano & Eyles 2003; Lavigne et al. 2010), and to predict the effects of climate change (e.g. Sulis et al. 2011; Levison et al. 2014a, 2014b).

The hydrogeological systems of eastern Canada have undergone major changes since the last deglaciation. During the accelerated isostatic rebound period that followed the Laurentide Ice Sheets retreat, new emerging recharge zones and increased potentiometric heads favoured a large invasion of meltwater into the newly formed Quaternary granular aquifers (Person et al. 2007). Major clay deposits from the Champlain Sea that followed the last deglaciation have isolated several areas of the lowland fractured bedrock aquifers. Since then, this ‘old’ water has been diluted at varying rates depending on recharge rates and the hydrogeological setting. In these conditions, the isotopic residence time measured at a given well will represent the travel time of the groundwater, but also the remaining signal of the initial paleo-recharge. Only limited knowledge on the ability of advective transport modelling to estimate groundwater travel times is available in these high recharge and high water table areas.

The objective of this work was to evaluate the benefits and limitations of using calculated isotope-derived groundwater travel times and major ion geochemistry to validate a calibrated regional groundwater flow model in the humid continental climate of southern Québec (Canada). A three-dimensional (3D) regional-scale steady-state groundwater flow model was built for the Centre-du-Québec region in southern Québec, using field data from aquifer characterization studies performed in the Bécancour River (Larocque et al. 2013), the Nicolet River and the lower Saint-François River watersheds (Larocque et al. 2015). The model was calibrated using heads and baseflows in MODFLOW (Harbaugh 2005). Available groundwater ³H/³He and ¹⁴C residence time data (Vautour et al. 2015; Saby et al. 2016) and water chemistry data (Meyzonnat et al. 2016) were then used to assess the quality of the calibration.

Study site

Hydrology and meteorology

The study area corresponds to the administrative region of Centre-du-Québec, which is located approximately midway between Montréal and Québec, on the southern bank of the St. Lawrence River (Figure 1). It covers an area of 7452 km² and includes the watersheds of the Bécancour and Nicolet rivers, as well as several small watersheds that discharge directly into the St. Lawrence River (Figure 1). The topography varies between the plain, with elevations of between 7 and 100 metres above sea level (m asl), and the Appalachian Mountains, with peaks of 700 m asl and several deep valleys with steep slopes.

Figure 1. Study area in the Centre-du-Québec region (Québec, Canada), with the locations of the two watersheds, of baseflow measuring stations, of wells sampled for ¹⁴C and ³H/³He (from Vautour et al. 2015 and Saby et al. 2016), and of wells sampled for major ions (from Meyzonnat et al. 2016 and Saby et al. 2016).

The climate of the region is characterized by cold winters and warm, humid summers. The long-term average air temperature varies from 6.2°C on the plain to 4.4°C in the Appalachians. Minimum air temperature (−15°C) occurs in January, while maximum air temperature (26°C) is in July. The average annual precipitation varies from 1018 mm in the lowlands to 1213 mm in the Appalachians. Precipitation occurs mostly as rainfall between April and November, while snowfall in the other months represents 23% of the total precipitation (Environment Canada 2016; station nos. 7028441, 7020305, 7022160 and 7025440).

Geology

The regional fractured aquifer is composed of rocks belonging to two geological provinces: the Appalachian Mountains in the southeastern portion of the basin, and the St. Lawrence Platform in the northwestern portion, both of Cambro–Ordovician age (Figure 2a). Geographically, these two geological provinces are part of the St. Lawrence Lowlands.

Figure 2. Geology of the study area: (a) bedrock geology (modified from Globensky 1987 and Slivitzky & St-Julien 1987), (b) bedrock confinement zones (modified from Larocque et al. 2013, 2015), and (c) geological cross-section (from Saby et al. 2016).

The fractured bedrock formations of Ordovician age are the main aquifer of the area. These can be separated into two major rock type zones: the sedimentary and the metasedimentary (Globensky 1987; Slivitzky & St-Julien 1987). The sedimentary rock zone covers the plain from the St. Lawrence River to the Appalachian Piedmont, and consists of limestone and carbonate-rich shales, shale, and mixed shale and fine sandstone. This area is slightly deformed. A network of southeast-dipping sub-horizontal faults is present. The metasedimentary rock zone mainly covers the Appalachian Mountains and consists of a wide range of schist, quartzite and phyllades. These rocks of the Appalachian Mountains are highly deformed by a network of faults and folds.

Unconsolidated Quaternary fluvioglacial deposits cover the fractured bedrock aquifer (Lamothe 1989), and can form superficial aquifers of limited aerial and vertical extent. Basal deposits are tills from the last two Quaternary glaciation episodes (Bécancour Till [early Wisconsinan] and Gentilly Till [Middle to late Wisconsinian]), followed by glacio-lacustrine, sandy and organic deposits. A thick clay layer deposited during the Champlain Sea episode (12–9.8 ka BP; Bolduc & Ross 2001) following the last deglaciation covers sandy deposits over a strip of 10 to 30 km parallel to the St. Lawrence River (Godbout et al. 2011; Lamothe & St-Jacques 2014). This clay creates confining conditions for the underlying fractured bedrock aquifer. Clay thickness is greater in the Nicolet River watershed (> 40 m) and tends to be lower in the northeast part of the study site (Bécancour watershed). The central part (between the clay plain and the Appalachian Piedmont) is composed of sand, patches of clay, and outcropping till and shale, creating a heterogeneous and unconfined to semi-confined hydrogeological system. In the Appalachian Mountains, reworked till and bedrock outcrops leave the fractured aquifer unconfined in its main recharge area (Figure 2b). In most of the valley bottom, silty material creates semi-confined to confined conditions.

A 3D geological model of the regional aquifer was constructed from borehole interpolation and a surface geological survey (Larocque et al. 2013, 2015). According to the model, overburden thickness varies between 90 m in the northeast part of the study site to less than 1 m in the Appalachian Mountains. Quaternary granular aquifers are not spatially extensive, and are generally located in complex stratigraphic systems consisting of sand, clay, till and varve sequences. Clay silt and till create semi-confined and confined conditions for the fractured bedrock aquifer.

Hydrogeology

The regional-scale fractured bedrock aquifer is used for drinking purposes by approximately half of the 115 municipalities present in the area. Larocque et al. (2013, 2015) compiled hydrogeological data from consultant reports and from the provincial wells drillers’ database (PWDD) managed by the provincial government (MDDELCC 2013). The piezometric map of the fractured bedrock aquifer (Figure 3a) shows a general groundwater flow pattern from the Appalachian Mountains to the St. Lawrence River. This regional flow is interrupted locally where the topographic gradient is steeper and where rivers deeply incise the surface (i.e. the downstream portions of the Bécancour and Nicolet rivers). The rivers are disconnected from the bedrock aquifer in several areas where thick clay deposits occur. The average water table depth is 4.8 m below the ground surface.

Figure 3. Hydrogeology of the study area: (a) piezometric map, and (b) distributed recharge and groundwater model limits (both panels modified from Larocque et al. 2013, 2015).

The Quaternary granular aquifers are discontinuous in the investigated area. Most of the municipal wells exploit unconfined surficial sand aquifers that are generally less than 10 m thick. In some areas, the granular aquifer can be up to 20 m thick, but its spatial extent is limited.

Available hydrogeological, isotopic and geochemical data

Reliable groundwater potentiometric data (from private consultant reports, as well as from data obtained by Larocque et al. 2013, 2015) are available for 535 wells (Figure 3a), and lower quality data (PWDD) are available for 11,240 wells. All the potentiometric data were used as calibration targets for the groundwater flow model.

Larocque et al. (2013, 2015) also conducted hydraulic tests to estimate aquifer transmissivity in 20 observation wells drilled into the bedrock aquifer (30 to 50 m of open boreholes; Figure 4). Results indicate that, despite substantial heterogeneity, there is no indication of a trend or a clear link between the aquifer lithology and the hydraulic properties. Results from the hydraulic tests show hydraulic conductivities ranging from 3.7 × 10⁻⁹ to 8.1 × 10⁻⁶ m.s⁻¹. Interpretation of specific capacity data from the well drillers’ database (Huntley et al. 1992; Richard et al. 2016) suggests a geometric mean hydraulic conductivity of 1.2 × 10⁻⁶ m.s⁻¹. Packer test results available for nine of the 20 observation wells indicate that horizontal hydraulic conductivity can vary vertically by two orders of magnitude (from 6.5 × 10⁻⁸ to 9.6 × 10⁻⁵ m.s⁻¹), but do not show any systematic trend with depth. Larocque et al. (2013, 2015) identified no hydrogeological influence of faults. Malgrange and Gleeson (2014) also showed the limited influence of the Appalachian fault system on the local hydrogeology.

Figure 4. Measured (interval) and mean (black dots) hydraulic conductivities (from Larocque et al. 2013, 2015): (a) for the fractured aquifer, and (b) for the granular aquifers.

Larocque et al. (2013, 2015) also conducted hydraulic tests in the granular deposits, and compiled results from consultant reports. Results from the hydraulic tests on 13 piezometers and 16 municipal wells show hydraulic conductivities ranging from 2.3 × 10⁻⁷ to 5.0 × 10⁻³ m.s⁻¹ (Figure 4b).

Recharge was estimated using a physically based and spatially distributed water budget model (Larocque et al. 2013, 2015). The model uses the Runoff Curve Number method (USDA 2004) to simulate runoff, while evapotranspiration is calculated using the Oudin et al. (2005) empirical equation. The model is calibrated to reproduce total runoff and baseflow simulated using the global MOHYSE hydrological model (Fortin & Turcotte 2007). Spatial information, such as soil type, land use and slope, was reported on a 500 m × 500 m grid. For each cell, a water balance was then compiled by successively subtracting runoff and evapotranspiration (ETP) from daily precipitation and snow melt. Because recharge estimation is based on a surface model and does not integrate stratigraphy, a correction factor is applied to semi-confined and confined areas to limit the recharge rates. For the semi-confined area, 40% of the potential recharge is assumed to reach the aquifer (Figure 3b). This is based on the hypodermic runoff/total flow ratio obtained from a surface flow model reported in Larocque et al. (2013, 2015), and is also in the range of that used by Rivard et al. (2014; between 9 and 52%), based on surface geology and confining conditions. No recharge was used in the areas where the aquifer is confined. The average recharge for the entire study area, calculated for the 1989–2009 period, is 156 mm.yr⁻¹ (Larocque et al. 2013, 2015).

Baseflow data were also used as calibration targets in the model. Baseflows were estimated by the provincial hydrological service using streamflow time series and the Eckhardt (2005) filter, as reported by C. Poirier (MDDELCC-Québec 2012). Flow rate time series are available for nine gauging stations (see Figure 1 for locations), and their duration varies between 20 and 80 years. Over these periods, baseflows represented between 24 and 39% of total flows at the stream gauging stations. Manual stream flow measurements were taken on seven small rivers (see Figure 1 for locations) using a Swoffer flowmeter (three repeated measurements for stations located on the Bécancour River watershed, and only one measurement for stations located on the Nicolet River watershed) during the low-flow periods of July and August between 2009 and 2012 (Larocque et al. 2013, 2015).

A regional groundwater chemistry survey was performed by Larocque et al. (2013, 2015), with further analysis and interpretation reported in Saby et al. (2016) for the Nicolet River watershed and Meyzonnat et al. (2016) for the Bécancour River watershed. Major and minor ion concentrations, as well as water chemistry types, are available for 211 wells distributed over the entire study area. Groundwater is dominated by the Ca-HCO₃ type in the Appalachian Mountains and in the unconfined areas of its piedmont. In the central and lower portion of the study area, cationic exchange between Ca and Na generates a Na-HCO₃ groundwater type, whereas in confined areas close to the St. Lawrence River, low recharge and diffusion with clay pore seawater from the Champlain Sea invasion (11 kyrs ago; Occhietti & Richard 2003) create an Na–Cl groundwater type.

Groundwater residence times were estimated using ³H/³He and ¹⁴C methods (Vautour et al. 2015; Saby et al. 2016). Groundwater samples were obtained from 36 private, municipal and observation wells during the summers of 2010 and 2013. Among these wells, three were installed in the granular aquifers and 33 in the fractured aquifer. Eighteen samples were analyzed only for ³H/³He, eight were analyzed only for ¹⁴C and 10 were analyzed for both. Water samples for ³H analysis were collected using 1 L Nalgene® bottles filled and sealed before shipment to the Environmental Isotope Laboratory (EIL) at the University of Waterloo. Liquid scintillation counting (LSC) is used by the EIL to quantify tritium. The detection limit is 0.8 TU. ³He was analyzed at the noble gas laboratories of the University of Michigan and the University of Tokyo using noble gas isotope mass spectrometry (see Vautour et al. 2015 for details; and also see Saby et al. 2016). ¹⁴C was analyzed at the Beta Analytic Laboratory in Florida (Nicolet River watershed samples) and at the Institut de Chimie des Substances Naturelles-Centre National de Recherche Scientifique in Gif-sur-Yvette, France (Bécancour River watershed samples).

The ³H/³He residence time was calculated taking care to separate the different sources of ³He in groundwater other than that produced from tritium decay (tritiogenic ³He or 3He_tri). These are: terrigenic helium (³He_terr) from mantle and/or crustal production; atmospheric helium dissolved at solubility equilibrium at the recharge (³He_eq); and atmospheric helium in excess of the solubility equilibrium amount (the so-called ‘excess air’ produced by air bubbles at the air/water table interface or ³He_ea; e.g. Schlosser et al. 1989). Data were reported on a Weise plot (Weise & Moser 1987), and equations from Schlosser et al. (1989) were used to determine the amount of ³He_tri. To test the validity of the obtained residence time, the initial ³H content prior to decay (i.e. ³H+³He_tri) was compared with the ³H in precipitation at the Ottawa Global Network of Isotopes in Precipitation (GNIP) station. Most samples fall within the calculated yearly-averaged Ottawa tritium input curve (see v 2016 and Vautour et al. 2015 for details on tritium dating).

Because of the occurrence of several pools of dead carbon in the studied aquifer, particularly in the Bécancour River watershed where dissolved methane occurs (Moritz et al. 2015), uncorrected and adjusted ¹⁴C groundwater residence times were obtained using the NETPATH-WIN^® software (El-Kadi et al. 2011), which is the only software suitable for methanogenic aquifers (Aravena et al. 1995). For wells in the Nicolet River and the lower Saint-François River watersheds, where methane was not directly measured, NETPATH-adjusted residence times were very similar to those from the classical Fontes and Garnier equilibrium model (Fontes 1992). This latter method was preferred by Saby et al. (2016) in this area. Only one well yielded reliable results using the Tamers model (Tamers 1975), probably because of its high levels of dead carbon (Saby et al. 2017).

In most cases, adjusted ¹⁴C residence times are very close to the uncorrected residence times. However, in a few cases they are almost half the uncorrected value. Because of the numerous dead carbon pools that are difficult to take into account (see Vautour et al. 2015 for details), it was decided to report both uncorrected and corrected residence times to give a range of possible ¹⁴C residence times in the study area (Figure 5).

Figure 5. Comparison of residence times estimated from ³H/³He and ¹⁴C with particle travel times. The ³H/³He residences times are compared with the minimum travel time, while the ¹⁴C-derived residence times are compared with the maximum travel time. The dashed circle represents wells in confined and semi-confined conditions. Both uncorrected and corrected (Netpath) ¹⁴C values are shown and are linked by the black line.

The mean ³H/³He groundwater residence time is 34 years (ranging from 3 ± 0.14 to 60 ± 1.4 years). The average uncorrected ¹⁴C groundwater residence time is 6234 years (varying between 226 ± 11 and 22 600 ± 1130 years). The average adjusted residence time is 4224 years (varying between 147 ± 7 and 17,050 ± 850 years).

Young ³H/³He groundwater residence times are observed in high-elevation recharge areas and in superficial granular aquifers. Older ¹⁴C groundwater residence times are found in confined and semi-confined areas in the lower part of the study site close to the St. Lawrence River and Lake Saint-Pierre (a widening of the St. Lawrence; see Figure 1). They are also observed in valley bottoms of the Nicolet River, where thick deposits of fine sediments create confining conditions (see Vautour et al. 2015 for details; and also see Saby et al. 2016).

Model description

Model geometry

The groundwater flow model was developed using the MODFLOW software package (Harbaugh 2005). Groundwater flow in the fractured bedrock and in the Quaternary granular aquifers was simulated in steady state, assuming a regional-scale flow behaviour similar to that of an equivalent porous medium. The simulated domain covers the entire study area and is composed of a 250 m × 250 m cell grid with 12 vertical layers (for a total of 186,560 cells). The elevation of the top layer corresponds to the surface topography, obtained from the provincial digital elevation model (10 m vertical resolution; MERN 2013). This first layer includes all Quaternary deposits, and extends to the top of the bedrock with a minimum thickness of 2 m and a maximum thickness of 80 m. The 11 other layers are parallel to the bedrock topography. The model base elevation is flat and has a constant slope, with an elevation ranging from −300 m in the lower part of the model (St. Lawrence River) to −120 m in the upper part (Appalachians). Layer thicknesses are 10 m (layers 2 to 4), 20 m (layer 5) and 40 m (layers 6 to 12). The average model thickness is 320 m. This arbitrary choice was a compromise between attaining reasonable model computation time and considering data availability. Little or no hydrogeological data were available below 200 m. Tran Ngoc et al. (2014) compiled data from the petroleum industry and show permeability data from 700 m deep close to 1 × 10⁻¹⁴ m²/s combined with the presence of brackish water. The base of the fresh groundwater system is thus considered to be located above 700 m.

The simulated area was separated into eight hydraulic conductivity (K) and effective porosity zones (same zonation for K and porosity). Model layer 1 was separated into three zones corresponding to the unconfined, semi-confined and confined areas (see Figure 2a). Layer 2 consists of one zone that represents the more fractured shallow bedrock. Layers 3 to 6 were separated into two zones delineated by the transition of sedimentary/low metamorphic to metamorphic rocks (Figure 2b). Finally, layers 7 to 12 were also separated into two zones, similarly to layers 3 to 6. K values were calibrated for the different zones. Effective porosity values were not calibrated, because almost no data are available for the study area. The effective porosities used in the model vary between 0.05 (for the shallow bedrock) and 0.015 for the deeper bedrock layers. These values are within the 1–5% interval identified by Tran Ngoc et al. (2014) for the Ordovician fractured aquifer. For the unconsolidated sediments, effective porosity values of 0.15, 0.08 and 0.01 were used for the unconfined (medium to fine sand), semi-confined (fine sand and silt) and confined (clay) zones respectively, similar to values used by Benoit et al. (2011).

Boundary conditions

The lateral boundaries of the study area (southeast, northeast and southwest) were set as no-flow limits. The downgradient limit that corresponds to Lake Saint-Pierre and the St. Lawrence River was represented with a constant head (elevation 4 m) and assigned to all the layers. Small lakes located in the studied watersheds were also simulated in layers 1 to 3 with a constant-head boundary condition at elevations corresponding to that of the DEM. The river network was generated with the ‘flow accumulation’ tool in ArcGIS (ESRI 2016), using the same 250 m × 250 m DEM as for topography. This was done to ensure that the river boundary conditions were correctly located in the topography. The resulting river network corresponds to a flow accumulation of 30 cells (i.e. the smallest watersheds cover an area of 1875 km²). The river system was represented in MODFLOW using a DRAIN boundary condition. This representation is well adapted to a shallow water table aquifer, where rivers and streams are generally fed by groundwater and where rivers rarely feed the aquifer (Gleeson et al. 2011). Drain elevations were set to that of the DEM. Conductance values were calculated from the ratio between the K value and the cell thickness, multiplied by the cell width (as suggested in MODFLOW). Thus, for a cell thickness of 5 m and a cell width of 250 m, the river conductance varied between 6 × 10⁻⁴ m².s⁻¹ (K = 1.2 × 10⁻⁵ m.s⁻¹) and 6 × 10⁻⁶ m².s⁻¹ (K = 1.2 × 10⁻⁷ m.s⁻¹). To avoid limiting flow to the drains, conductance values were set to 1.2 × 10⁻³ m².s⁻¹ for all the drain cells. The spatially distributed 20-year average recharge (1990–2010) from Larocque (2013, 2015) was applied to the top layer. An evapotranspiration boundary was applied to the first layer of the model to limit the water table from being above the soil surface as there is no evidence of widespread artesian areas in the study site. An evapotranspiration rate of 500 mm.y⁻¹ and an extinction depth of 1 m were used. This boundary condition also helps to deal with the possible overestimation of recharge, and provides information on groundwater seepage areas.

Model computing, model calibration and particle tracking

The model was computed using the NWT Newton Package (Niswonger et al. 2011) of MODFLOW-2005. Inter-cell conductance was solved using the Upstream-Weighting (UPW) package. This package treats the non-linearity of cells drying and rewetting using a continuous function for groundwater head, helping to lessen convergence problems.

The groundwater flow model was calibrated using a trial-and-error procedure, by manually adjusting the hydraulic conductivity in the seven K-zones. A calibration was also done with PEST (Doherty 2016). The head calibration targets included all available head measurements (from consultant reports, from values reported in Larocque et al. [2013, 2015], and from the well drillers’ database [MDELCC 2013]). Calibration targets also included the baseflow data extracted from flow rate time series, and from field measurements. Effective porosities and the spatially distributed recharge were not calibrated.

Groundwater travel times were calculated using MODPATH (Pollock 2016). One hundred particles were assigned to each of the 36 wells for which isotopic residence times were available. The particles were equally distributed among layers 2 through 6 (10 to 70 m below the surface) according to the open borehole length of each well. Particles were backtracked to their source and forced to pass through a weak sink (a particle stops in a cell if a sink represents more than 99% of the water balance of the cell). The same method was used for the 211 wells with water type data. Average borehole depth (40 m) was used and particles were assigned to layers 2 to 4. The total travel time of each particle was used to calculate statistics and verify the presence of statistically significant differences between water types using the JMP software (SAS Institute Inc. 2012).

Results

Model calibration and groundwater flow conditions

The manual calibration showed that all the simulated heads were equally distributed along the 1:1 line, not showing any systematic bias related to elevation (Figure 6). The mean error (ME) on heads from consultant reports and from Larocque et al. (2013, 2015) data was 1.7 m, while the mean absolute error (MAE) was 4.1 m. The root mean square error (RMSE) was 5.9 m and the normalized RMSE was 1.4% (i.e. less than the 10% threshold generally recognized as indicative of a reliable calibration). These statistical results indicate a good fit between observed and simulated heads. The PEST calibration (scattergram not illustrated) showed little improvement over the manual calibration, with an ME of 2.0 m and an MAE of 3.8 m. The RMSE was 5.4 m and the normalized RMSE was 1.3%.

Figure 6. Measured and simulated heads obtained via manual calibration. Mean error (ME), mean absolute error (MAE), root mean square error (RMSE), and normalized RMSE (NRMSE) are also shown. Data from MDDELCC (2013) are shown, but statistics are calculated using only data from Larocque et al. (2013, 2015).

The baseflows simulated by the manually calibrated model were lower than the average measured values for most of the stations, except in the cases of four manual measurement stations (Figure 7). In general, the baseflows estimated from available time series were better simulated than those from manual measurements, with an MAE of 150,667 m³.d^-−1(the absolute error was based on the mean value). Simulated baseflows from the PEST simulation show the same pattern as for the manual calibration, with a higher MAE, of 152,383 m³.d⁻¹.

Figure 7. Measured or estimated baseflows compared to simulated baseflows. The intervals represent minimum and maximum values from baseflow separation or from multiple field-measured values and the dots represent the mean value.

The calibrated hydraulic conductivities of the sediments (layer 1) were within the range of measured values (cf. Figure 4a), from 5.8 × 10⁻⁵ m.s⁻¹ for sand or unconfined areas, 1.16 × 10⁻⁶ m.s⁻¹ for silt and sand or semi-confined areas, and 1.16 × 10⁻⁸ m.s⁻¹ for clay and silt or confined areas. The calibrated hydraulic conductivity of the fractured shallow bedrock (layer 2) was 5 × 10⁻⁶ m.s⁻¹. K values were calibrated to be 1.2 × 10⁻⁷ m.s⁻¹ for layers 3 to 6 (same K value for sedimentary/low metamorphic and metamorphic zones; see Figure 2a), and were thus within the range of the available field-measured values (cf. Figure 4b). K values were calibrated to be 5.8 × 10⁻⁸ m.s⁻¹ for layers 7 to 12 (same K value for sedimentary/low metamorphic and metamorphic zones (see Figure 2a). No measurements were available for the deeper bedrock layers. The calibrated values suggest a decrease in hydraulic conductivity with depth. A similar decrease in K values with depth was obtained by Lavigne et al. (2010), while a more rapid decrease was reported by Sanford (2017). A single vertical anisotropy value of 1 (K_h/K_v) was calibrated manually for all zones, except for the confined zone for layer 1 (calibrated anisotropy of 10). The calibrated K values from the PEST runs were very similar to those from the manual calibration for layers 1 to 6. They were lower (1.2 × 10⁻⁹ m.s⁻¹ for the metamorphic zone and 6.6 × 10⁻⁹ m.s⁻¹ for the sedimentary/low metamorphic zone)) for layers 7 to 12. The PEST-calibrated vertical anisotropy remained the same for layers 1 to 6 (Kh/kV = 1), and increased to 100 in the two zones of layers 7 to 12.

The evapotranspiration boundary condition removed 52% of the total recharge. Simulated heads were close to the surface in these areas, creating seepage zones and leading to substantial uptake by the evapotranspiration boundary condition. The resulting net simulated recharge (i.e. recharge minus flux to the evapotranspiration boundary) over the entire modeled area was 121 mm.yr⁻¹, compared with 156 mm.yr⁻¹ reported by Larocque et al. (2013, 2015). The areas of overestimation (i.e. where the evapotranspiration boundary condition was active) generally corresponded to the valley bottoms of the Appalachian Mountains. The proportion of groundwater discharge to the stream is 47% of the total recharge. The PEST-calibrated model results in a higher uptake by the evapotranspiration boundary (105 mm.yr⁻¹) and a lower net recharge (106 mm.yr⁻¹). The proportion of groundwater discharge to the stream is reduced to 41% of the total recharge.

Particle tracking results

For each of the 36 wells for which groundwater residence times were estimated based on isotopic tracers, the model provided a probability density function of travel times from the array of particles reaching the well (Figure 8). The frequency distribution of simulated travel times is highly variable, ranging from 2 to 18,174 years. Most distributions have an inflection point at a cumulative frequency of approximately 0.2, which represents the transition between the second layer, with higher hydraulic conductivity (fractured bedrock), and the lower layers (3 to 6), with lower hydraulic conductivity.

Figure 8. Cumulative frequency distribution of travel times for the 36 wells for which isotope-derived residence times were available. The dotted line indicates the frequency distribution for all wells.

For most of the wells for which estimated ³H/³He groundwater residence times were available (Figure 9), the simulated travel times were within the range of the dating method (max ≈ 60 years). The manually calibrated model overestimated the ³H/³He groundwater residence times by 256 years on average, while the PEST calibration resulted in an overestimation of 160 years on average. At a cluster of four points (dashed ellipse on Figure 5), the simulated travel times were greater than 100 years. Among these wells, three were in the confined bedrock aquifer in the lower portion of the Bécancour River and the Nicolet River watersheds, and one was in the Piedmont of the Bécancour River watershed, where the bedrock aquifer is semi-confined.

Figure 9. Comparison of the simulated groundwater travel times with the dominant water types found in the study area. Median, mean, 25th and 75th percentiles are presented.

The simulated travel times generally underestimated measured values for wells with uncorrected measured ¹⁴C residence times, by 4401 years on average. The maximum simulated travel time was 14,934 years (i.e. similar to the maximum isotopic residence time of 22,600 ± 1130 years), occurring in an area where the bedrock aquifer is confined, close to Lake Saint-Pierre. The largest underestimations (by 7618, 8955 and 9781 years) of ¹⁴C residence times were for three wells located in the valley bottom of the Appalachian Mountains. These wells also contained ³H/³He, which indicates recently recharged groundwater with isotope-derived groundwater residence times of 48, 32 and 56 years, respectively. The simulated travel times also generally underestimated corrected ¹⁴C residence times, with an average error of 2073 years. The correction improves the results for most of the wells. As for the uncorrected values, the largest underestimations (6069, 7618 and 9791 years) occurred for wells located in valley bottoms of the Appalachian Mountains. The PEST-calibrated model generally underestimated uncorrected measured ¹⁴C residence times, by 1210 years on average, and overestimated corrected ¹⁴C residence times, with an average error of 1661 years.

Comparing simulated mean travel time and groundwater types revealed that these are a good indicator of the overall aquifer dynamics (Figure 9). Ca-HCO₃ groundwater (associated with recharge areas) was associated with shorter travel times on average than Na-HCO₃ or Na–Cl groundwater water types (with longer water–rock interactions). The median travel times were 68, 274 and 738 years for the Ca-HCO₃, Na-HCO₃ and NaCl water types, respectively. The differences in travel times between the three water types were statistically significant at α = 0.01 (Student t-test, P value < 0.0001). The PEST calibration also results in a statistically significant difference, at α = 0.01 (Student t-test, P value < 0.0001), between mean travel times of the three water types, but with higher median residence times. PEST calibration results in median travel times of 244, 843 and 1617 years for the Ca-HCO₃, Na-HCO₃ and NaCl water types, respectively (not shown in Figure 9).

Discussion

Regional groundwater flow

The results of the comparison between observed and simulated groundwater levels and between observed and simulated baseflows indicated that, overall, the model was able to represent the steady-state hydrodynamics of the study area. Error statistics were comparable to those from other groundwater flow modelling studies (Castro & Goblet 2003; Lavigne et al. 2010; Levison et al. 2014a; Rivard et al. 2014), and similar for the manual calibration and the PEST calibration. Errors on simulated heads could be caused by local-scale heterogeneities.

The discrepancies between manually measured baseflows and simulated baseflows could result from the fact that manually measured flows and the available flow rate time series were not representative of the 20-year average recharge conditions. More specifically, measurements of the Nicolet River watershed were made during a prolonged dry period, which may explain model overestimation (white dots without error bars in Figure 7). Underestimation of the baseflows for the manual measurements on the Bécancour River watershed (white dots with error bars in Figure7) can be explained by simplification of the stratigraphy. The aquifers associated with these manually measured stations are mainly overlain by confining units throughout the watershed in the model, whereas in reality a surficial aquifer is present and contributes to maintaining baseflows. A similar effect of stratigraphic simplification on the simulated baseflows was also observed by Juckem et al. (2006) in a study of the Driftless Area of Wisconsin. They showed that baseflow in a small watershed (< 50 km²) was more sensitive (i.e. increased discrepancy) to hydrostratigraphic simplification than in larger watersheds. The low recharge rates applied to confined areas induce low simulated baseflow. Increasing river and drain conductance did not improve results. This is acceptable, however, because baseflow time series obtained using filtered total flow rates generally overestimate actual groundwater discharge. Croteau et al. (2010) and Rivard et al. (2014) obtained similar underestimation of measured baseflows with a steady-state groundwater flow model. An overestimation of baseflows occurred in unconfined areas, while an underestimation occurred in the clay plain. The underestimation of the baseflow could also be caused by the spatial resolution used to represent the river network, or the grid resolution near the drain boundary conditions used to represent rivers and streams (Haitjema et al. 2001; Brunner et al. 2010). It is important to underline that, similarly to the heads, the baseflows were simulated well with the manual calibration and the PEST calibration.

A large flux of water is removed from the model by the evapotranspiration boundary conditions, especially in valley bottoms where no drains were present. The level of detail of the drain network probably influenced this volume of water, because several small valley bottoms are not represented as having drains. The evapotranspiration boundary removed groundwater that in reality is routed to the stream network. A smaller amount of groundwater-simulated recharge compared to that estimated with surface infiltration modelling was found in this study, similar to the results of Rivard et al. (2014). It is important to note that a large part of the difference corresponds to water that infiltrated at high elevations, and was removed through the evapotranspiration boundary in the valley bottoms. Because recharge overestimation occurred only in these locations, the spatially distributed recharge used in other areas of the model can be considered a satisfactory estimate of the actual average recharge rate.

In terms of the water budget, discharge to streams and rivers represents approximately half (47%) of the total water balance. This substantial removal of infiltrated water by the river and stream network can have large implications in terms of groundwater vulnerability and contamination. For example, short groundwater flowpaths between infiltration and discharge in a river imply that natural attenuation to reduce concentrations of contaminants sourced at the surface might not have the time to occur. Short flowpaths also indicate that contaminant dilution might not be a significant process in reducing concentrations. Recent measurements of contaminant concentrations (nitrate, pesticides and pharmaceutical compounds) in groundwater in the study area might reflect this groundwater flow dynamic (Saby et al. 2017).

The discharge of groundwater to Lake Saint-Pierre and the St. Lawrence River was very low, representing only 0.01% (0.02 mm.yr⁻¹) and 0.17% (0.28 mm.yr⁻¹) of the water budget, respectively. These small fluxes show that regional groundwater flow to the St. Lawrence is almost nonexistent. Groundwater fluxes between the two main watersheds were equivalent to a very small net flux, of 0.59 mm.yr⁻¹, from the Nicolet River watershed to the Bécancour River watershed. This indicates that although the hydrological limits do not correspond exactly to the hydrogeological limits, the two boundaries are very similar.

Groundwater travel times

The results obtained here show a relatively good simulation of ³H/³He-based groundwater residence times and an underestimation of some of the ¹⁴C-based residence times. This could suggest that the model does not provide an adequate portrait of the variety of groundwater travel times. Sanford et al. (2004) observed similar patterns when comparing ¹⁴C and particle travel times in a regional groundwater flow model, and suggested that the discrepancy with variation in past recharge might explain these results. The underestimation of ¹⁴C residence times in recharge areas was explained by slower-than-calibrated recharge in the recent past, and the overestimation of ¹⁴C residence times in discharge areas was explained by higher-than-calibrated recharge in the more distant past. However, this suggests the presence of an important regional flow component, which the current results do not support in the study area. Another explanation could be the presence of groundwater originating from the large recharge that followed the deglaciation 12 kyrs ago (e.g. Person et al. 2007). The mixing of old groundwater, infiltrated at that time, with more recently recharged water could explain the discrepancy between ³H/³He- and ¹⁴C-derived residence times. This process was also suggested by Vautour et al. (2015), Saby et al. (2016) and Méjean et al. (2017) to explain the occurrence of mixed residence times in the study area, where old groundwater with low ¹⁴C activity and high radiogenic ⁴He content mixed with more recent recharge, enriched in tritium and with atmospheric helium.

Sanford (1997) and Wassenaar and Hendry (2000) suggested that ¹⁴C stored in paleo-pore water could diffuse into the faster flowpath, leading to an underestimation of ¹⁴C-derived travel time by advective transport. This mechanism cannot be ruled out in the study area, although only one sample (with a high ¹⁴C residence time, of 17 kyrs) had a clear Na-Cl chemistry, suggesting that ¹⁴C had exchanged with pore seawater from the Champlain Sea clays. In the study area, the major problem in correcting ¹⁴C residence times was the occurrence of the additional dead carbon reservoir of methane (Moritz et al. 2015; Vautour et al. 2015). The methane content was measured in 130 wells in the study area as part of a shale gas study (Pinti et al. 2013), but only a few correspond to those presented here. It cannot be ruled out that ¹⁴C residence times could be younger for several of the wells reported in Figure 5 due to this additional methane-¹⁴C.

Other aspects of the conceptual model could explain the discrepancy between measured and simulated travel times. As shown in the results section, wells that have groundwater travel times higher than those calculated by the ³H/³He method are located in semi-confined or confined aquifers. The stratigraphic simplification in these zones can induce artificially low recharge in areas where punctuated recharge exists in otherwise confined units. In such cases, including dispersion would not improve the ability to model short travel times, except if dispersion is sufficiently large to make the particle reach the closest recharge zone. This could be the case for wells located in semi-confined aquifers, because these areas are more discontinuous and frequently alternate with unconfined zones.

At the regional scale, the variation in effective porosity within each hydraulic conductivity zone is not considered to have a large impact on simulated groundwater travel times compared to the spatial distribution of recharge within the same unit (Portniaguine & Solomon 1998). Effective porosities are rarely measured for fractured bedrock aquifers. Improving the quality of simulated travel times would necessitate a finer spatial discretization of K values than the zonation used in the model. This would be difficult to justify given the available field-measured data. With an average discrepancy of 4401 years (2073 years for corrected residence times) between isotopic ¹⁴C residence times and the particle tracking travel times, the effective porosity would need to be increased by more than an order of magnitude to improve results. The porosity values would thus be greater than the range expected from the work of Tran Ngoc et al. (2014), of 0.5–10%. However, as suggested by Sanford (1997), a dual porosity system could greatly impact ¹⁴C residence time estimation by allowing the diffusion of older ¹⁴C water from the low-flow layers to the high-flow layers. This could lead to younger ¹⁴C residence times and a better fit with advective travel times.

In Figure 5, the divergence from the 1:1 line is related to recharge and confining conditions at each well, and to the actual mixing between the ³H/³He, younger water and ¹⁴C, older water end members. Saby et al. (2016) showed that groundwater in the study area results from the mixing of a younger component, containing post-bomb ³H and ¹⁴C, and a pre-bomb fossil water component, containing natural background ³H and dead ¹⁴C. The amount of the pre-bomb component in the mixture was found to be as high as 98% in some areas (Saby et al. 2016). These particular mixing conditions could explain the overestimation of ³H/³He-obtained residence times by the model in highly confined areas (dotted ellipse in Figure 5), because the regional simplification may have removed small recharge areas. These results also suggest that a recalibration of the model to ¹⁴C residence times by reducing hydraulic conductivity or recharge rates would have led to an overestimation of the travel time in recharge areas (i.e. an overestimation of ³H/³He residence times). Including the presence of so-called ‘old water’ in a steady-state regional groundwater flow model represents a real difficulty, and the comparison of particle travel times with ¹⁴C-derived residence times should be made with caution and interpreted only qualitatively. Szabo et al. (1996) suggested that calculated groundwater residence times should be considered qualitatively rather than used directly as calibration targets for simulated travel times, thus serving instead as a validation of the conceptual model.

Although the manually calibrated model provided similar errors on heads and baseflows to the PEST-calibrated model, it was clearly superior in simulating adequate mean travel times for the three groundwater types. The PEST-calibrated model simulated unrealistically long travel times for the Ca-HCO₃, Na-HCO₃ and NaCl water types. This could be due to the higher vertical anisotropy (Kh/kV = 100, for both K zones) and the lower K of the deep bedrock (layer 7–12, for both K zones) resulting from the PEST automatic calibration.

Model sensitivity

The sensitivity analysis was first based on the results from PEST to identify the parameters that have the largest influence on travel times. The most influential parameters were then modified manually to compare their effect on heads’ MAE and on mean travel times. The overall parameter sensitivity determined from PEST indicates that the hydraulic conductivity of the fractured bedrock (first 10 m, layer 2) is the parameter to which the model is the most sensitive, followed by that of the shallow bedrock (layers 3–6) of the St. Lawrence platform (Figure 10a). Sensitivity of K the fractured bedrock (layer 2) reflects the connectivity of the aquifer to the river network and the importance of baseflows in the model. The model is also sensitive to the sediments’ K. Vertical anisotropy of the shallow bedrock (layers 3–6) in the metamorphic zone is also among the most sensitive parameters.

Figure 10. (a) parameter sensitivity from the PEST analysis. Numbers in brackets indicate the model layers, and (b) variations in the mean groundwater travel times of the three groundwater types and in the MAE on heads resulting from changes in recharge, hydraulic conductivities, and vertical anisotropies. In panel (b) the codes used on the x-axis are explained in the legend below the figure.

In light of the PEST-identified model sensitivities, the six most sensitive parameters were selected and were modified manually to estimate their specific impacts on heads and on travel times. Impact of variations in recharge rate was also tested. The results indicate little impact of parameter changes on head calibration statistics (Figure 10b). However, when looking at the impact of parameter changes on travel times (Figure 10b), it is clear that using head calibration statistics is insufficient to evaluate the quality of the model calibration in terms of groundwater types. Varying hydraulic conductivities for fractured bedrock (layer 1) or shallow bedrock (layer 2–6) of the St. Lawrence platform resulted in an unrealistic range of travel times for each of the groundwater types. Low hydraulic conductivity values increased the mean travel time for Ca-HCO₃ groundwater to greater than 1000 years, which is conceptually unrealistic for groundwater associated with recharge areas. A similar observation was made for travel times of Na-Cl groundwater (less than 1000 years, after increasing hydraulic conductivities for layers 3–6 of the sedimentary/low metamorphic zone). Sanford (2011) reported similar challenges when simulating advection travel times in synthetic models. This author also showed that travel times are sensitive to the spatial variation in recharge. This may explain the low sensitivity of travel times to the variation in recharge, as the tests conducted here were only done on the recharge rate and not on its spatial distribution. Because porosity does not influence the water budget in steady-state simulations, the impact of varying porosity on the travel times is expected to be proportional to the change in porosity values. Because this impact is straightforward, it was not tested here, but would be expected to be much lower than that of the other tested parameters.

Conclusions

This study demonstrated the difficulty of modelling groundwater residence time with particle tracking in a context of high water table and important mixing between groundwater recharged after the last deglaciation and recently recharged groundwater. The results indicated that young isotopic water residence times are more easily simulated than old groundwater, because of the highly dynamic groundwater system, which is closely linked to the river network through baseflow. In parallel, a good simulation of old ¹⁴C residence times, representing groundwater infiltrated several thousand years ago (and not linked to a long travel time), is mostly achieved in confined areas that are located far from the recharge zones, where no mixing with newly recharged water occurs.

Although the PEST-calibrated model provided similar errors on heads and baseflows to the manually calibrated model, and a small improvement in groundwater travel times, it generates older travel times for three groundwater types. The validation of the calibrated model with residence time tracers and particle tracking thus appears to be essential.

Results from this study indicate that groundwater types are useful to validate groundwater flow models. Short travel times should be associated with recharge-type groundwater (Ca-HCO₃) and longer travel times with confined-type groundwater (Na–Cl), and travel times of the Na-HCO₃ water type should lie between these two. This combination of groundwater type based on major ion chemistry and particle travel times has rarely been used to date in groundwater modelling.

The regional geological simplifications used in this work are realistic and do not create a bias in the model. However, at the local scale, this simplification can generate problems, such as the underestimation of baseflows for small watersheds. The use of a spatially distributed recharge obtained from a surface water budget plays an important role in buffering the impact of the geological simplification. This modelling work has allowed the integration of data from regional-scale groundwater characterization projects. Such integration permits the validation of individual analyses from different research domains, such as geochemistry, isotopic geochemistry and hydrogeology. Despite the difficulty related to the mixing of different groundwater masses, the use of isotopic tracers and groundwater types allowed the validation of a regional groundwater model. These tracers should not be directly included as calibration targets, but should instead be used in model validation. Because major ion groundwater types are relatively easy to measure, and low cost, they should be included in groundwater model validation. Further work is now needed to better understand the transient effect of recharge on groundwater travel times at the regional scale, and the modelling of paleo-recharged groundwater in regional models.

Acknowledgments

The authors would like to thank the Quebec Ministry of Environment (Ministère du Développement durable, de l’Environnement et de la Lutte contre les changements climatiques) and the partners of the regional groundwater assessment projects of the Bécancour and the Nicolet-Saint-François area for their financial contributions. We would also like to thank all the people who were involved in the two regional groundwater projects (research associates, MSc students, field assistant and technicians). Finally, we wish to thank the reviewers for their comments that greatly improved the quality of this work.