Regional groundwater flow dynamics and residence times in Chaudière-Appalaches, Québec, Canada: Insights from numerical simulations

Abstract

As part of the Quebec PACES III provincial groundwater resources assessment programme (Programme d’acquisition des connaissances en eaux souterraines), a regional-scale two-dimensional numerical groundwater model was developed in the Chaudière-Appalaches region, Québec, Canada. The model considers groundwater flow, transport of groundwater age and the influence of a fault on the flow system and its implications for groundwater quality. By including deep and shallow flow systems, the study helps fill a knowledge gap with respect to intermediate flow systems and the role they would play during potential energy resource development including shale gas exploitation from the Utica Shale. Physical and chemical hydrogeological data, including an analysis of ¹⁴C in dissolved inorganic carbon in sampled groundwater, supported a regional conceptual flow model forming the basis for numerical simulations. The numerical model is first calibrated to regional piezometry through a semi-automated workflow using the inverse model PEST. Although some evidence for deeper regional flow exists, the area appears to be dominated by local flow systems on maximum length scales of about 5 km, with significant flow through the top 40 to 60 m of the fractured sedimentary rock aquifer. This regional-scale flow model is also supported by the local hydrogeochemical signatures. Simulated mean groundwater ages show young shallow water of < 100 years with rapid increases in age with depth suggesting diffusion-controlled age evolution. Groundwater age is likely being perturbed in the vicinity of the Jacques Cartier River fault, which can act as both a barrier and a preferential pathway, provided permeability contrasts with the surrounding rock are at least two orders of magnitude.

Dans le cadre du projet Programme d’acquisition de connaissances sur les eaux souterraines III pour la région de Chaudière-Appalaches, Québec, au Canada, l’étude présente une analyse approfondie de l’influence des dynamiques d’écoulement sur la qualité des eaux souterraines dans un contexte régional. L’écoulement régional, le transport d’âge et l’impact d’une faille sur la qualité de l’eau souterraine sont étudiés par l’entremise de modèles numériques bidimensionels. En considérant les systèmes d’écoulement profonds et peu profonds, l’étude clarifie les dynamiques intermédiaires qui les relient et le rôle qu’elles pourraient jouer dans l’exploitation potentielle des gaz de shale de l’Utica. La combinaison de connaissances hydrogéologiques physiques et chimiques, y compris une analyse des concentrations de ¹⁴C dans les eaux souterraines échantillonnées, a conduit à l’ébauche d’un modèle conceptuel de l’écoulement régional sur lequel des simulations numériques ont été basées. Le modèle est d’abord calibré à l’aide du logiciel PEST en comparant les charges simulées à la piézométrie régionale. Bien que le modèle numérique affiche l’existence d’un écoulement régional profond, la région à l’étude apparaît être dominée par des systèmes d’écoulements locaux à des échelles maximales d’environ 5 km, avec un écoulement significatif dans le roc fracturé peu profond. L’écoulement actif se limitant à une profondeur maximale de 40 à 60 m du roc fracturé. La simulation du transport advectif-dispersif de l’âge montre des eaux jeunes de < 100 ans près de la surface et une augmentation rapide de l’âge avec la profondeur suggérant l’importance de la diffusion comme processus de transport. Enfin, les résultats suggèrent que, si le contraste de perméabilité entre la faille et la roche mère est au moins de deux ordres de grandeur, la faille de la Rivière Jacques Cartier pourrait agir à la fois comme une barrière et une zone d’écoulement préférentielle.

Introduction

Understanding relationships between local, intermediate and regional groundwater flow systems (as defined by Tóth 2009) is critically important for management of water resources as well as to reduce the risk of contamination from deep energy resource development (e.g. shale gas or geothermal energy). In this context, a regional-scale groundwater research programme was initiated in the Chaudière-Appalaches (CA) region, Québec, as part of the groundwater resources assessment programme (Programme d’acquisition des connaissances sur les eaux souterraines, PACES) of the Quebec Ministry of Environment (Ministère du Développement durable, de l’Environnement et de la Lutte contre les changements climatiques; MDDELCC 2017). The present study was part of one of many PACES projects, which focused on regional groundwater flow dynamics and natural hydro-geochemistry and geochemical processes in the Chaudière-Appalaches region (Figure 1; referred to here as PACES-CA).

Figure 1. Location and topography of the Chaudière-Appalaches study area, and location of the two-dimensional regional cross-section model, modified from Lefebvre et al. (2015).

The CA study area, south of Quebec City, spans over 14,625 km² and has a population of approximately 276,000, mainly concentrated near major waterways including the St. Lawrence River, and the Chaudière River and its tributaries (Lefebvre et al. 2015). Although surface water is abundant, almost half of the water usage (42%) in the CA region is extracted from the subsurface.

As groundwater is fundamental to the local economy, groundwater resources within the CA region, in terms of both quantity and quality, have been the focus of a number of studies over the past three decades, culminating in the current PACES project. These projects have helped to characterise the resource, focusing on the most densely populated Chaudière River watershed (McCormack 1982; Rousseau et al. 2004; Benoît et al. 2011; Brun Koné 2013; Benoît et al. 2014; Carrier et al. 2014; Lefebvre et al. 2015). Based primarily on data collected from shallow residential wells (with average depths of approximately 60 m), these studies have provided important information on groundwater quality and aquifer properties, concluding that groundwater flow occurs primarily within permeable valley sediments and in the top fractured part of the bedrock (<30 m; Brun Koné 2013; Benoît et al. 2014; Crow and Ladevèze 2015; Ladevèze et al. 2016; Ladevèze 2017). Furthermore, Benoît et al. (2014) and Lefebvre et al. (2015) found geochemical signatures left behind by the invasion of the Champlain Sea, which flooded the northern areas of the CA region (Figure 1) about 12,000 to 10,000 years ago (Parent and Occhietti 1988).

In parallel, the potential exploitation of unconventional gas resources from the Utica Shale of the St. Lawrence Platform geological province, which includes parts of the CA region, led to a series of additional geological and hydrogeological studies led by the Geological Survey of Canada and the INRS (Institut national de la recherche scientifique). These studies aimed to assess the risks of groundwater contamination associated with shale gas exploitation through upward fluid migration and to establish the baseline hydrogeological conditions prior to exploitation. Focusing their case study on the region of Saint-Édouard-de-Lotbinière (Figure 1), Bordeleau et al. (2015) and Rivard et al. (2014) found, despite local flow regimes being more apparent, indications of relatively deep regional flow emerging near the Jacques-Cartier River fault line (Figure 2). Their findings support the conceptualisation of regional groundwater flow as a Tóthian flow system (Tóth 1999, 2009), in which local flow regimes are embedded within larger intermediate or regional flow regimes (Goderniaux et al. 2013). As shale gas exploration areas border the St. Lawrence River (Lavoie et al. 2014; Rivard et al. 2014; Moritz et al. 2015), which is a clear outlet for regional groundwater flow, broadening the understanding of how such systems might interact with local flow systems becomes even more important.

Figure 2. Conceptual model for the Chaudière-Appalaches regional groundwater model, showing (a) flow and geochemical evolution with water types G1, G2 and G3; and (b) the geologic cross-section of the deep subsurface (after Séjourné et al. 2013).

The aim of the broader PACES-CA study was to reinforce understanding of the regional flow system in the Chaudière-Appalaches region, and included investigating the characteristic time scales of groundwater flow processes through a review of the regional geochemistry and from estimates of groundwater age based on isotopic evidence (Lefebvre et al. 2015). The current paper focuses on testing a conceptual model of the regional hydrogeological system using a two-dimensional numerical flow and age transport model. The numerical modelling aims to quantify the magnitude of regional flow, estimate the maximum depth of active groundwater flow, investigate the distribution of groundwater residence times, and determine the influence of major normal faults on regional flow.

Methodology

The study objectives are met by first carrying out a thorough review of the regional geological context and Quaternary history along with a synthesis of the state of knowledge of the regional geochemistry and an analysis of groundwater ages. Combining the physical and chemical hydrogeological knowledge led to the development of a regional conceptual flow model.

The conceptual model was then tested by simulating the regional flow system within a two-dimensional vertical plane, about 70 km long, oriented roughly southeast–northwest from a regional topographic flow divide in the Appalachians towards the St. Lawrence River. All flow simulations are run using the FLONET/TR2 numerical code (Molson and Frind 2017), assuming steady-state, saturated isothermal conditions. The flow model is first calibrated using the PEST inverse calibration model (version 13.6; Doherty 2015). The calibrated flow model and various fault scenarios are then applied in the TR2 transport model to simulate the distribution of groundwater ages. The role of faults is investigated by testing various fault configurations based on the findings of a hydromechanical characterisation of the Saint-Édouard-de-Lotbinière area conducted by Ladevèze (2017) and Bordeleau et al. (2015, 2018).

Study area

Geologic and geographic context

Regional landscape and topography are characterised by a combination of various geological settings, dominated by the geological provinces of the Appalachians and the St. Lawrence Platform (Figure 2b). Although both geologic provinces are primarily sedimentary in origin, the Appalachians are highly deformed and slightly metamorphosed whereas the St. Lawrence Platform has undergone almost no deformation.

Situated towards the south of the study area, the Appalachians are the remains of a weathered mountain range and are characterised by a hilly terrain. The Appalachian region is marked by granitic intrusions forming small mountains, including the Petits Monts Mégantic and the Monts Notre Dame mountains (Figure 1) reaching elevations up to 850 and 1000 m, respectively (Slivitzky and St-Julien 1987). The intrusions are surrounded by a belt of small hills, which are bordered to the south by the Appalachian plateau and to the north by a foothill area that marks the junction with the sedimentary plains of the St. Lawrence Lowlands (Slivitzky and St-Julien 1987). With a history of high tectonic activity, the Appalachian region is marked by the presence of lineaments oriented NE–SW which manifest primarily in the form of major faults (Caron 2012; Lefebvre et al. 2015).

In the northernmost part of the study area, the St. Lawrence Platform consists of a sedimentary succession sitting unconformably on the Canadian Shield (Globensky 1987). Due to tectonic compression, normal faults mark the gradual lowering of the sedimentary sequence, which plunges below the Appalachian Basin (Castonguay et al. 2010) as shown in Figure 2b.

The Cambrian–Ordovician sedimentary sequence of the St. Lawrence Platform is a typical marine transgression–regression sequence resulting from the opening and closing of the Iapetus Ocean. From the bottom up, the sequence begins with the Potsdam Group, a poorly cemented sandstone formation becoming well cemented in its upper portion (Cloutier et al. 2006). Above this formation sits the commonly named carbonate platform, which is separated into four groups including dolomitic sandstone, dolostone, limestones and shales (Cloutier et al. 2006; Séjourné et al. 2013).

At depths ranging from 800 to 3000 m, the Utica Shale, which can be considered marl due to its high carbonate content, has a potential for shale gas production that was tested by vertical and horizontal exploration wells, some of which were subjected to hydraulic fracturing (Thériault 2012; Séjourné et al. 2013; Lavoie et al. 2014; Bordeleau et al. 2015). The overlying Lorraine group, composed of silty shales with mostly non-calcareous sediments interbedded with fine sandstone and clayey siltstones (Thériault 2012) is the top formation of the St. Lawrence Platform sedimentary succession over most of the study area.

Many past studies have identified the regional aquifer as being located in the uppermost portion of the bedrock where fracture density is highest (Benoît et al. 2014; Lefebvre et al. 2015; Ladevèze et al. 2016). Major faults could also represent areas where the level of fracturing could be relatively high, although the role of these faults in the regional flow system is not clear.

Surficial sediments throughout the area are on average 3–5 m deep but can reach 30 m in the valleys. As the study area is generally covered with till deposits of various thicknesses and assemblies, the nature of the till units defines the degree of confinement of the bedrock aquifer. In the Appalachian Highlands, surficial sediments are thin or non-existent, therefore unconfined conditions prevail for most of the fractured bedrock aquifer, with the exception of valleys where the aquifer is generally semi-confined (Lefebvre et al. 2015). Moreover, glaciation–deglaciation processes, the presence of the Champlain Sea approximately 10,000 years ago, and the formation of post-glacial lakes and shorelines led to heterogeneous deposits with varying proportions of sand and silt. These sediments together with alluvial deposits represent significant infiltration pathways when directly in contact with the bedrock aquifer. On the other hand, scattered silt and clay sediments, deposited at the bottom of stagnant water bodies following glacial retreat, prevent recharge and are responsible for the accumulation of organic material (Lefebvre et al. 2015). These low-permeability deposits are present mostly within the area formerly covered by the Champlain Sea, specifically in the northwestern area of the St. Lawrence Lowlands.

Total precipitation over the study area is estimated to be 1149 mm/year, with an average of 166 mm/year of recharge reaching the fractured rock aquifer (Lefebvre et al. 2015). Spatial distribution of recharge is highly variable. In the St. Lawrence Lowlands, areas covered with thick layers of fine sediments can receive less than 15 mm/year whereas rock outcrop areas can receive more than 300 mm/year. In the Appalachian Highlands, where the surficial sediments are fairly thin and permeable, aquifer recharge tends to be higher, varying from 100 mm/year to more than 300 mm/year. Additional studies on specific watersheds within the study area have reported overall recharge rates to the bedrock aquifer of between 27 and 38 mm/year (Benoît et al. 2014 and L’Heureux 2005, respectively). Discrepancies between recharge values could be attributed to its spatial variability.

Regional-scale conceptual model

Combining the geologic interpretation with the regional piezometric and hydrogeochemical data presented in Lefebvre et al. (2015) led to development of an initial regional conceptual flow and geochemical evolution model (Figure 2a). Most of the groundwater sample data reflected a rapid transit time, thus confirming that flow occurs mostly at shallow depths in the top portion of the fractured rock aquifer where it is under open system conditions with respect to the atmosphere.

Based on the analysis of groundwater samples taken from residential wells (Janos 2017), three major water types were found to prevail in the study area, named G1, G2 and G3. First, G1 samples, with typical recharge water geochemical signatures (Ca-HCO₃) are found in both the St. Lawrence Lowlands and the Appalachian Highlands (Figure 2a). Most of the samples of this water type were identified as being in an open system with respect to the atmosphere, while those in a closed system have the lowest radiocarbon ages found among the three water types. G1 waters would correspond to the active flow layer.

In the St. Lawrence Lowlands, the presence of G2 (Na-HCO₃)-type waters indicates that salty water which had infiltrated the aquifer during the Champlain Sea invasion has been replaced by fresh water. Above a mean groundwater age threshold of 10,000 years, the salty water would have been replaced by fresh recharge at least once. Therefore, the waters sampled above this threshold should have a geochemical signature corresponding to G1 or, if cation exchange sites with Na+ ions are present, to G2 water which corresponds to the first step of the freshening process.

G3 water types (Na-HCO₃/Na–Cl) are related to stagnation zones, long travel paths or mixing with older waters such as formation brines. These water samples have a maximum radiocarbon age of 11,680 years BP, and support the hypothesis that some wells pump from older water horizons. Bordeleau et al. (2015, 2018) collected two water samples in the Saint-Édouard-de-Lotbinière area, near the northern end of the study area, at a depth of 48 m below the surface without disturbing the water column in the well, and dated the samples to over 1.5 Ma (million years ago). The work of Vautour et al. (2015) and Saby et al. (2016) also emphasises the presence of old recharge water and connate brines in near-surface groundwater samples collected in the St. Lawrence Platform.

In the St. Lawrence Lowlands, the horizontal hydraulic gradients are lower than in the highlands, thus leading to slower groundwater flow rates, which is reflected in water samples showing signs of water freshening. The presence of such waters in shallow residential wells indicates that the interface between the active flow zone and slow-moving groundwater is close to the surface. Furthermore, the presence of waters having old radiocarbon ages and long flow paths, and possibly mixing with formation brines from the Appalachians, points to a possible regional flow path emerging near the St. Lawrence River. Overall, the geochemical and isotopic evidence (Lefebvre et al. 2015) suggests that the flow system in the CA region seems to follow the Thóthian regional groundwater flow concept, which is characterised by nested local, intermediate and regional flow systems (Tóth 2009).

Additional insight into developing the two-dimensional (2D) conceptual model for this study was provided by a three-dimensional (3D) numerical flow model of the Chaudière River watershed (Brun Koné 2013; Benoît et al. 2015). Using a maximum model depth of 300 m, including shallow unconsolidated sediments and the fractured sedimentary rock, this 3D model showed that most groundwater flow is focussed in the shallow rock. However, the model did not investigate deep regional flow. Additional groundwater flow models in two and three dimensions focusing on flow systems through shallow fractured rock have also been developed for other regions in Québec (Basses-Laurentides: Nastev et al. 2005; Châteauguay: Lavigne et al. 2010; Montérégie Est: Laurencelle et al. 2013; Saint-Charles River wastershed: Cochand 2014; Graf 2015; Outaouais: Montcoudiol et al. 2017). Although important insight was gained through these studies, the importance and role of deep regional flow in these flow systems remains to be defined.

Numerical modelling

Objectives

The numerical model presented herein is specifically aimed at better understanding the relationships between local, intermediate and regional flow systems within the CA region. Furthermore, this study aims at identifying how the main geologic features may influence the behaviour of regional flow. In particular, the study area is marked by the presence of major faults with uncertain hydraulic characteristics and which are presumed to be located within a regional flow discharge zone. A better understanding of the role of faults in regional and deep groundwater flow is particularly important regarding the risk to shallow groundwater quality from potential development of unconventional hydrocarbons (Rivard et al. 2014; Lefebvre 2017).

Numerical modelling is used here to investigate the role of these faults in groundwater discharge. While the role of a conduit-type fault has been investigated for a generic regional basin with a geologic context comparable to the St. Lawrence Platform (Gassiat et al. 2013), the hydraulic behaviour of fault zones in the Chaudière study area has yet to be characterised (some early work was presented by Ladevèze et al. 2016; Ladevèze, 2017).

Modelling strategy

In order to best capture the dynamics of regional flow in the CA region, a 2D groundwater flow model was developed following a regional flow path (Figure 1), which was then coupled to an advective-dispersive age transport model. The flow and age transport simulations were completed using the finite element FLONET/TR2 code (Molson and Frind 2017).

In FLONET/TR2, the 2D steady-state groundwater flow equations are defined as a function of hydraulic heads (Equation 1) and stream functions (Equation 2) (Molson and Frind 2017):(1) (2)

where x and y are the horizontal and vertical coordinate directions, respectively (L), K_xx and K_yy are the principal components of the hydraulic conductivity tensor (L/T), is the hydraulic head (L) and ψ is the stream function (L²/T).

The steady-state saturated flow model was calibrated using PEST version 13.6 (Doherty 2015), a model-independent parameter estimation code. Following calibration, a full sensitivity analysis of selected parameters was completed (not shown here; for details see Janos 2017). The calibrated flow model was then used to simulate the transport of age mass. Finally, a parametric study on the role of faults was performed in order to investigate their potential influence on groundwater flowpaths.

Domain description

The numerical model domain corresponds to a 2D projection of the geologic cross-section and conceptual groundwater flow model. The model follows an inferred regional flow path which extends from a local peak elevation in the Appalachian Highlands to the St. Lawrence River, being 69 km long with a maximum depth of 9.1 km. Three distinct geological provinces are included in the model: the Grenville crystalline bedrock, relatively flat-lying units of the sedimentary St. Lawrence Platform, and the Appalachians. The entire domain is covered by a discontinuous layer of surficial deposits of marine and glacial or fluvioglacial origin.

Steady-state regional flow model

Assumptions and limitations

Several assumptions have been made in order to simplify the CA model and reduce over-parameterisation. First, given the regional scale of the study, fractured rock is represented as an equivalent porous medium. Also, using a 2D model assumes that groundwater flow is perfectly aligned with the cross-section and thus transverse hydraulic and age gradients are neglected. Since the cross-section is parallel to regional topographic and piezometric gradients, it is considered well representative of the regional flow direction.

Furthermore, the flow system is simulated under steady-state conditions and thus does not account for seasonal (or climatic) recharge variations. Transient recharge is expected to have the most impact on local flow systems, with little or no significant impact on regional flow. Furthermore, water table fluctuations induced by seasonal recharge variations are also expected to be negligible with respect to inherent uncertainty, including the errors associated with the extrapolation of the regional potentiometric map, which is assumed to represent an annual average water level. Variations over geological time (e.g. ice advances and retreats) are also neglected.

Moreover, the model assumes that groundwater in the domain has a uniform temperature, viscosity and density. The average thermal gradient in the area is between 15 and 20°C/km (Grasby et al. 2011) which, based on an estimate of the Rayleigh number for average conditions found in the study area, is not sufficient for the onset of thermal convection (Hiscock and Bense 2014). Also, although measured total dissolved solids concentrations in the study area reach up to 16,785 mg/L (Bordeleau et al. 2015, 2018) and up to 340,000 mg/L in the deep formations of the St. Lawrence Platform (Tran Ngoc et al. 2014), the effects of density and viscosity gradients on flow are assumed to be less important compared to the 5–6 order of magnitude logarithmic decrease with depth of the hydraulic conductivities (discussed later).

Domain extent and discretisation

The initial top boundary elevations of the model, representing the water table, were extracted from the regional potentiometric map (Lefebvre et al. 2015), while the elevations of the bedrock interface and the nine distinct surficial deposit layers were obtained from a detailed cross-section of surficial deposits traced roughly along the modelled section (Lefebvre et al. 2015). For modelling purposes, the rock formations within the cross-section were divided into five different zones representing different rock types: the St. Lawrence Platform, the Appalachian Highlands, the Utica Shale, the carbonate platform and silicate bedrock of the Grenville Province.

The finite element flow model is discretised using linear triangular elements with a horizontal discretisation of 50 m and a vertical discretisation between 0.1 and 50 m (increasing from top to bottom), for a total of 358,000 nodes and 797,640 elements.

Material properties

Material properties represented by the model can be classified as surficial deposits, the fractured bedrock aquifer and deep units.

Surficial deposits

Surficial deposits along the cross-section are composed of nine different hydrostratigraphic units (Lefebvre et al. 2015): organic deposits, alluvium, lakeshore delta sediments, glaciomarine delta deposits, coarse and fine-grained littoral glaciomarine deposits, fluvio-glacial deposits, till and undifferentiated sand. Each unit is represented by a layer of deformed elements in the model (Table 1). All layers of the surficial deposits were assumed to have an average porosity of 0.3.

Table 1. Physical description of the surficial deposits, including reported and calibrated hydraulic conductivities.

Mod. unit	Geologic group/description	Hydraulic description	Reported log K (m/s)			Reference	Calibrated log K (m/s)
org	Organic deposits (phagnum and ericaceous peat)	Poorly drained		−4.52		Tecsult (2008)	−4.73
				−5.00		Brun Koné (2013)
				−6.48		Ladevèze (2017)
all	Alluvium (silt and sand)	Permeable		−4.00		Tecsult (2008)	−6.00
			−6.00	to	−3.00	Brun Koné (2013)
			−7.44	to	−4.06	Ladevèze (2017)
			−9.00	to	−4.70	Domenico and Schwartz (1998)
lgd	Lakeshore delta (sand, silts, gravel)	Permeable	−6.00	to	−3.00	Brun Koné (2013)	−6.00
			−6.70	to	−3.52	Domenico and Schwartz (1998)
mgd	Glaciomarine delta (sand, silts, gravel)	Permeable	−6.00	to	−3.00	Brun Koné (2013)	−5.9
			−6.70	to	−3.52	Domenico and Schwartz (1998)
mgbg	Littoral glaciomarine (coarse) (gravel, sand, silts)	Permeable		−4.00		Tecsult (2008)	−4.4
			−6.00	to	−3.00	Brun Koné (2013)
			−8.12	to	−5.27	Ladevèze (2017)
			−6.70	to	−3.52	Domenico and Schwartz (1998)
mgbf	Littoral glaciomarine (fine) (silts, sand, gravel)	Relatively permeable		−4.00		Tecsult (2008)	−5.97
			−6.00	to	−3.00	Brun Koné (2013)
			−8.12	to	−5.27	Ladevèze (2017)
			−6.70	to	−3.52	Domenico and Schwartz (1998)
grav	Glaical/fluvioglacial (sand and gravel)	Very permeable		−3.52		Tecsult (2008)	−3.96
			−6.00	to	−3.00	Brun Koné (2013)
			−6.70	to	−3.52	Domenico and Schwartz (1998)
till	Glacial tills (sand, silt, clay)	Spatially variable		−8.52		Tecsult (2008)	−6.05
				−8.70		Tecsult (2008)
				−6.22		Brun Koné (2013)
			−12.00	to	−5.70	Domenico and Schwartz (1998)
sab	Sand	Permeable		−6.52		Tecsult (2008)	−4.86
				−7.00		Brun Koné (2013)
			−9.00	to	−5.00	Ladevèze (2017)
			−8.10	to	−3.52	Domenico and Schwartz (1998)

Till is the most prominent unit over the study area and it is thus expected that its hydraulic conductivity will control groundwater recharge into the fractured rock aquifer. The till unit also has the most heterogeneous composition; typically very sandy in the area of the cross-section, its texture can vary abruptly from a permeable sand facies with a few pebbles to an essentially impermeable stony and clayey silt (G. Légaré Couture, personal communication 2013). Uncertainty with respect to the hydraulic conductivity of the till layer was addressed through calibration and sensitivity analysis (Janos 2017).

Fractured bedrock

Hydraulic conductivities of the bedrock were obtained through the updated PACES-CA database (Lefebvre et al. 2015) from which the extensive public well drillers’ log data were used to estimate hydraulic conductivities (Figure 3). These data originate from specific capacity data while accounting for drilling and testing biases (M. Laurencelle, personal communication, June–July 2015; Lefebvre et al. 2015). Calculated hydraulic conductivities in the shallow bedrock (within the top 100 m of the rock surface) are between 3.8 × 10⁻¹⁰ and 6.5 × 10⁻⁴ m/s in the St. Lawrence Platform (241 data points) and between 5.2 × 10⁻¹² and 3.8 × 10⁻³ m/s in the Appalachian Highlands (4966 data points). Although the hydraulic conductivity of the bedrock is highly variable for a given depth, spanning more than six orders of magnitude, a clear trend of decreasing hydraulic conductivity with depth can be observed for both geologic provinces (Figure 3; Lefebvre et al. 2015). The calibrated hydraulic conductivities were obtained by adjusting the parameters of an exponential decay function given by Equation (3):

Figure 3. Vertical distribution of measured and model-generated hydraulic conductivities in the Chaudière-Appalaches region. Grey dots represent estimated values from the well drillers’ log (Lefebvre et al. 2015). Black squares show the median measured value per 5-m intervals, and the horizontal whiskers show the 25th and 75th percentiles. Solid lines show the calibrated distribution of hydraulic conductivities in the fractured bedrock, and dashed lines show different conductivity scenarios tested in the sensitivity analysis of Janos (2017).

(3)

where z(m) is the depth below the top of the bedrock surface, K(z) is the hydraulic conductivity at depth z, K_min is the minimum hydraulic conductivity asymptotically approached by the curve, K_max is the maximum hydraulic conductivity, and α is a curve decay constant (values given in Figure 3).

Deep layers

The Grenville basement rock, the St. Lawrence carbonate platform and the Utica Shale are described in the literature as having distinct hydraulic properties (Table 2). In the model, these units are differentiated from the top fractured unit of the St. Lawrence Platform (Lorraine formation) and the Appalachian Highland formations. Conductivities assigned to these deep formations are given in Table 2.

Table 2. Description of the modelled geologic rock units, including the reported and modelled hydraulic conductivities and porosities (includes properties reported by Gassiat et al. 2013). Bedrock anisotropy obtained from Freeze and Cherry (1979). Bold values are used in the model.

Model unit	Geologic group	Literature description	log K (m/s) min to max	Por. (%)	Reference
Above shale in SLP	Lorraine	Fine-grained sedimentary^a	−10.21 to −6.71		Gleeson et al. (2011)
	Silty shales with mostly non-calcareous sediments interbedded with fine sandstone and clayey siltstones	Lorraine unit matrix Bécancour (eastern Canada)^a^,^b	−10.61	3.9	Tran Ngoc et al. (2014)
		Fractured rock^c	−6.47 to −4.96		Ladevèze (2017)
		Fractured rock in St. Lawrence Platform formation	−9.43 to −3.19		Lefebvre et al. (2015)
		Minimum values in model (log Kx/log Kz)	−9.01/−10.01	5
Above shale in AH	Humber zone	Fractured rock	−8.96 to −5.39		Ladevèze (2017)
	Carbonates, sandstones and shales	Fractured rock in Humber formation	−11.28 to −2.03		Lefebvre et al. (2015)
		Minimum values in model (log Kx/log Kz)	−8.98/−9.68	5
Utica Shale	Utica Shale	General^a	−15.21 to −9.21		Neuzil (1994)
	Calcareous black shale	General^a	−12.21 to −8.21		Freeze and Cherry (1979)
		General^a	−15.21 to −12.21		Flewelling and Sharma (2014)
		Utica (eastern Canada)^a	−14.26 to −8.31		Séjourné et al. (2013)
		Marcellus (Northeastern USA)	−7.91		Soeder (1988)
		Barnett	−9.21		Montgomery et al. (2005)
		Utica unit Bécancour matrix (eastern Canada)^a^,^b	−10.73	4	Tran Ngoc et al. (2014)
		log Kx/log Kz in model	−9.5/−11.5	5
Carbonate Platform	Trenton	Trenton unit Bécancour matrix (eastern Canada)^a^,^b	−7.81	3.4	Tran Ngoc et al. (2014)
	Clayey limestone interbedded with shales	Trenton unit Bécancour global (eastern Canada)^a^,^b	−5.10		Tran Ngoc et al. (2014)
	Black River	Limestone interbedded with sandstone
	Chazy	Clayey and sandy limestone
	Beakmantown	Beauharnois unit Bécancour matrix (eastern Canada)^a^,^b	−8.26	1	Tran Ngoc et al. (2014)
	Dolomitic sandstone becoming pure dolostone then dolomitic limestone	Beauharnois unit Bécancour global (eastern Canada)^a^,^b	−8.43		Tran Ngoc et al. (2014)
		Theresa unit Bécancour matrix (eastern Canada)^a^,^b	−8.43	2.6	Tran Ngoc et al. (2014)
		Theresa unit Bécancour global (eastern Canada)^a^,^b	−8.10		Tran Ngoc et al. (2014)
	Potsdam	Cairnside unit Bécancour matrix (eastern Canada)^a^,^b	−8.13	3.3	Tran Ngoc et al. (2014)
	Poorly cemented sandstone formation becoming well cemented	Cairnside unit Bécancour global (eastern Canada)^a^,^b	−6.59		Tran Ngoc et al. (2014)
		Covey Hill unit Bécancour matrix (eastern Canada)^a^,^b	−7.83	6.3	Tran Ngoc et al. (2014)
		Values in model	−8/−9	1
Grenville	Grenville	Crystalline^a	−9.21 to −4.71		Gleeson et al. (2011)
	Unfractured metamorphic and igneous	Unfractured metamorphic and igneous rock	−12.21 to −8.71		Freeze and Cherry (1979)
		Calibrated hydraulic conductivities in Outaouais	−12 to −6.	1–3	Montcoudiol et al. (2017)
		log Kx/log Kz in model	−11/−11	0.05

^a Calculated from permeability values using density = 1000 kg/m³, viscosity = 8.90 x 10⁻⁴ Pa·s and g = 9.81 m/s².

^b Median value given in study.

Flow boundary conditions

Vertical boundaries, corresponding to the St. Lawrence River to the northwest (left boundary) and to a regional topographic and piezometric peak to the southeast (right boundary), are assumed to represent symmetric physical boundaries of the flow system and were thus assigned a no-flow boundary condition (Figure 4). It is also assumed that no significant flow occurs below a depth of 9100 m; thus, a zero-flux condition was applied to the base of the model. At stream and river locations, the top boundary is constrained by imposed heads. Given that recharge is highly variable over the study area, the top boundary was separated into 295 sections, each with a unique recharge imposed as a Type-2 flux boundary condition which was determined with PEST during calibration. At these recharge locations, the water table elevations were determined iteratively through vertical grid deformation (maximum of 10 layers) until the computed water table head equalled its elevation within a tolerance of 10⁻⁴ m (Molson and Frind 2017).

Figure 4. Geological structure, hydraulic conductivities and boundary conditions for the numerical flow model.

Flow model calibration

The numerical flow model was calibrated through a semi-automated workflow using PEST (version 13.6, Doherty 2015). PEST is a model-independent parameter estimation code which uses the Gauss–Levenberg–Marquardt search algorithm (GLMA) to minimise the difference between a set of observed (known) values and the corresponding simulated results (Doherty 2015). This software can be particularly useful when a large number of parameters need to be adjusted (calibrated) and when the relationship between the parameters and the observations is not linear (e.g. Hayley et al. 2014; Siade et al. 2015; Zhu et al. 2015). Nevertheless, when field observations are not sufficient to appropriately constrain a model, calibration solutions will be non-unique and human input is needed to identify the most plausible solution. To this effect, PEST was used in this project to create many realisations of the flow model, each constraining the parameters of the model differently through various iterations of the parameter estimation. The realisation with the best fit to observations and with the most plausible parameters was selected as the base case calibrated model.

Calibration was completed by coupling PEST to FLONET and allowing PEST (1) to modify the hydraulic conductivities of each surficial deposit layer, (2) to define the decrease in hydraulic conductivity in the fractured bedrock using Equation (3), and (3) to vary the imposed Type-2 Neumann surface recharge boundary conditions. Given the spatial variability in recharge calculated from one-dimensional (1D) infiltration models (Benoît et al. 2014; Lefebvre et al. 2015), recharge values imposed along the top boundary of the model were allowed to vary between 0 and 900 mm/year.

The observations, to which PEST was trying to find the best fit, were taken to be the elevation of the water table at each of the top boundary nodes, and were extrapolated from a regional potentiometric map based on 15,790 observation wells throughout the basin (Lefebvre et al. 2015). Because these data were collected over several decades and at different times of the year, they are considered an accurate long-term representation of the water table and are consistent with the steady-state assumption in the model. The average recharge values imposed on the model by PEST were verified to be consistent with the values calculated from the 1D infiltration model, HELP (Benoît et al. 2014; Lefebvre et al. 2015), based on the approach of Croteau et al. (2010).

With respect to the observed heads, which in the study area vary from 0 to over 600 m, the final calibrated model heads have an absolute mean error of 2.06 m, a root mean square error (RMSE) of 3.28 m and a maximum absolute error of 18.26 m (Figure 5). The error on simulated heads is well below 5% (30 m) of the overall variations in heads in the model domain, which according to Anderson et al. (2015) defines the acceptability of a groundwater model hydraulic head calibration. A longitudinal profile comparing the observed and simulated heads, together with surface topography and recharge, is provided in Figure 6.

Figure 5. Model calibration: Observed heads (from interpolated regional piezometry) against simulated heads.

Figure 6. Calibrated recharge fluxes, simulated hydraulic heads, observed water table elevations (from interpolated piezometry) and topographic elevation along the model cross-section. Negative fluxes are exiting the model domain and positive fluxes are entering (corresponding to recharge).

The average overall calibrated recharge to the bedrock is within 1.04% of the average recharge estimated with the HELP model (Table 3; Figure 6). Despite some regional differences between the model-calibrated average recharge to bedrock and the average values calculated with HELP, these calibrated recharge rates remain within the range of values estimated with HELP (Table 3). The model can therefore be considered a fair representation of the natural groundwater system at the regional scale.

Table 3. Average recharge to the bedrock calibrated from the HELP model and from the current two-dimensional flow model (SLL = St. Lawrence Lowlands, AHL=Appalachian Highlands).

Average recharge to the bedrock	Modelled with HELP (Lefebvre et al. 2015)	Calibrated flow model	Absolute error (%)
Overall (mm/year)	166	164	1.04
SLL (mm/year)	85 (< 15 to > 300)	63	25.64
AHL (mm/year)	187 (100 to > 300)	326	74.27

Calibrated hydraulic conductivities of the surficial deposits are within limits found in the literature, with the exception of the till layer having a calibrated hydraulic conductivity slightly higher than reported (Table 1). The same trend is observed in the calibrated hydraulic conductivity distribution of the fractured bedrock (Figure 3), which is an order of magnitude higher than the median hydraulic conductivities inferred from the well drillers’ log data.

This discrepancy between measured and calibrated hydraulic conductivities can be attributed to the well-documented scaling effect (Schulze-Makuch et al. 1999; Nastev et al. 2004; Zhang et al. 2006). The smaller the measurement volume of a hydraulic test performed in an aquifer, the less likely it is to intercept a preferential flow pathway such as a high-K channel or interconnected fractures. Schulze-Makuch et al. (1999) observed that upscaling the radius of influence of hydraulic testing in a given heterogeneous aquifer could lead to an increase in measured hydraulic conductivity of up to one order of magnitude. This scale effect was also observed for fractured rock aquifers of the St. Lawrence Platform (Nastev et al. 2004).

In this study, hydraulic conductivities of the fractured bedrock (Figure 3) were mostly obtained from specific capacity tests for residential wells, which have a relatively small radius of influence of approximately 20 m (Nastev et al. 2004). Lefebvre et al. (2015) describe heterogeneities within the tills of the CA region which would lead to similar scale effects of hydraulic conductivity measurements.

Flow model results

Figures 7 and 8 show the simulated potentiometric field and streamlines (flow lines) of the calibrated groundwater flow model, which shows evidence of a Tóthian flow pattern with embedded regional, intermediate and local flow systems (Tóth 1999, 2009).

Figure 7. Calibrated steady-state flow model showing hydraulic heads and streamline distribution.

Figure 8. Calibrated flow model: (a) Zoom 1 (Appalachian Highlands), (b) Zoom 2 (St. Lawrence Lowlands close to St. Lawrence River), (c) Zoom 3 (St. Lawrence Lowlands close to Appalachian front), showing hydraulic heads and streamline distribution (location shown in Figure 7).

The simulated flow system is consistent with previous studies (Brun Koné 2013; Benoît et al. 2014) which suggested a dominance of local groundwater flow systems, mostly concentrated in the permeable sediments and within the top 20 to 40 m of the fractured bedrock aquifer. Although the simulated local flow systems in the Appalachian Highlands extend in some areas down to 800 m below the bedrock surface (Figures 7 and 8a), the relative significance of the volume of flow decreases rapidly with depth. These local flow systems have a maximum horizontal scale of approximately 5 km.

In this study, active flow is defined as being in those areas which receive an average groundwater flux of more than 1 mm/year, which is reached at a depth of approximately 30 m below the rock surface in the St. Lawrence Lowlands and approximately 60 m below the rock surface in the Appalachian Highlands (Janos 2017). The difference in the depth of the active flow zone between these two geological provinces highlights the effects of topography and recharge on the shape of the flow systems.

Furthermore, the prominent topographic features of the Appalachian Highlands also generate an intermediate flow system originating from the highest topographic peaks to the south and emerging at the Appalachian foothills, spanning more than 20 km (Figures 7 and 8a).

Under the assumed simulated conditions, the model also suggests the presence of a deep regional groundwater flow path which originates from the Appalachian Highlands and emerges near the St. Lawrence River (Figure 7). However, this regional system (defined between the two deepest streamlines in Figure 7), containing only about 0.5 m³/year per metre of transverse width, represents less than 0.005% of the total flow in the domain. Mean groundwater age simulations presented later also show that this deep flow would be extremely slow, with mean ages on the scale of geological time (tens of millions of years).

A sensitivity analysis (not shown; see Janos 2017) showed that the regional flow system is relatively insensitive to variations in the hydraulic conductivity of the surficial deposits, while in the Appalachian Highlands the depth of active flow seemed very sensitive to the hydraulic conductivity of the bedrock and to recharge variations. Deep intermediate and regional flow was highly affected by the bedrock hydraulic conductivity and to a lesser degree by the imposed recharge.

Groundwater age simulations

Modelling approach

In addition to using groundwater flownets, simulated mean groundwater ages can also be useful for interpreting flow systems. Groundwater age is defined as the time elapsed since water entered the groundwater system (Goode 1996; Kazemi et al. 2006).

In the TR2 groundwater model, mean groundwater age is taken to be the average (mean) age of a packet of water molecules at a given point in space (Bethke and Johnson 2008; Molson and Frind 2017), and is calculated using Goode’s (1996) method for direct simulation of groundwater age. The approach is more realistic than kinetic or advective age, since it allows age mixing (dispersion) by applying the advection-dispersion transport equation (Equation 4; see also Molson and Frind 2012):(4)

where A is the mean age (T), D_ij is the hydrodynamic dispersion tensor (L²/T), v is the groundwater velocity (L/T) and the term + 1 represents the age growth factor (1 day/day). Janos (2017) compares these advective-dispersive ages with particle-track based advective ages, and with measured ¹⁴C ages.

Transport model boundary conditions and parameters

The calibrated steady-state groundwater flow model was used as the base velocity field to simulate mean groundwater ages within the 2D cross-section of the CA regional model. A finite element mesh identical to the flow model was used in the age transport model.

In the age transport model, lateral and bottom boundaries are set as zero age-mass gradient Type-2 (Neumann) boundaries. Those sections of the top water table surface which act as inflow recharge boundaries are set as a Type-1 (Dirichlet) boundary condition with a fixed age of A = 0 days. Surface nodes acting as discharge points, i.e. nodes connected to elements with an outward-pointing velocity vector, are applied a zero-gradient age boundary condition. A growth rate of +1 day/day ensures the aging of water as the simulation progresses over time until steady-state is reached between the constant aging of water and rejuvenation from younger recharge water.

The dispersion tensor D_ij in Equation (4) is defined using the Lichtner formulation (Lichtner et al. 2002) which has four dispersivity components: two longitudinal (α_LH and α_LV) and two transverse dispersivities (α_TH and α_TV). Based on scaling relationships of Xu and Eckstein (1995) and Schulze-Makuch (2005), the longitudinal horizontal dispersivity of the current transport model was set to 50 m. The longitudinal vertical dispersivity (α_LV) was set to 10 m, while the two transverse dispersivities were set to 0.05 m. The molecular diffusion coefficient was 1 × 10⁻¹⁰ m²/s, similar to a conservative tracer. Similar dispersion and diffusion terms were applied in 2D and 3D age transport modelling by Montcoudiol et al. (2017) and Molson and Frind (2012).

The Peclet accuracy criterion, which constrains the spatial discretisation by limiting the ratio of advection to dispersion (Molson and Frind 2017), is respected throughout most of the domain. Using adaptive time-stepping, the Courant stability criterion (Molson and Frind 2017) is also generally respected for all time steps. Since the steady-state age simulations showed good convergence and relatively minor numerical oscillations, the spatial and temporal discretisation was considered acceptable.

Mean groundwater age distribution

The advective-dispersive groundwater age simulation reached steady state at 47.5 Ma for the base case scenario (Figures 9 and 10), which corresponds to the maximum simulated groundwater residence times in the Grenville basement rock. The nested flow systems create a unique mean age distribution which is organised into a complex pattern.

Figure 9. Simulated steady-state mean groundwater ages with superimposed streamlines from the calibrated flow model. (Plot enlargements for Zooms 1, 2 and 4 are provided in Figure 10; Zooms 3 and 5 can be found in Janos 2017).

The active flow layer within the local flow systems is characterised by short residence times, with mean groundwater ages being generally less than 100 years (Figures 9 and 10). Interestingly, the age simulation shows that the depth of the local flow systems extends well below the active flow layer. The maximum mean age of groundwater in these local flow systems is about 100,000 years (Figure 10a–c). While numerous ground surface discharge zones maintain relatively young ages within these shallow flow systems, mean ages nevertheless increase exponentially with depth which is attributed to the exponential decrease in the bedrock hydraulic conductivity, and to the natural distribution of regional flow, as shown by Vogel (1967) for simple single-compartment homogeneous systems with a single discharge zone. Goderniaux et al. (2013) presented similar behaviour in numerical simulations of nested flow systems.

Figure 10. Simulated mean groundwater ages with superimposed streamlines from the calibrated flow model: (a) Zoom 1, (b) Zoom 2, (c) Zoom 4. (Locations shown in Figure 9).

Furthermore, while the presence of an intermediate flow system in the Appalachian Highlands was clear from the flow simulation, the simulated mean groundwater ages also suggested the existence of such a system with a similar age distribution in the St. Lawrence Lowlands, between about 15 and 30 km along the section (Figure 9).

However, contrary to the Appalachian Highlands, in the St. Lawrence Lowlands the age distribution is marked by an abrupt transition between the regional and the smaller scale flow systems. In some portions of the cross-section, for example, the mean age gradient spans two orders of magnitude over very short vertical spatial scales on the order of 50 to 100 m (Figure 10b and ^{Figure 10}c). Numerous sharp vertical ‘plumes’ of older water also appear along the profile where deep flow is discharging to surface water (Figure 10). As it reaches the high-permeability shallow aquifers, this older water is rapidly diluted by faster flowing young water. This behaviour is consistent with Zijl’s (1999) observation that in real steady-state flow systems, mixing zones between transport systems are relatively thin, and thus the interface is generally sharp between water bodies with different water qualities (or ages) which belong to different flow systems. Where the age interface is horizontal, vertical diffusive age mixing prevails.

On a regional scale, the average residence times in the deep flow system gradually increase towards the discharge area (Figure 9). Regional flow in the St. Lawrence Lowlands seems to follow a deep preferential pathway through the carbonate platform before emerging near the St. Lawrence River. This flow pattern and age distribution is due to the difference in the assigned hydraulic conductivities between the Lorraine formation, the carbonate platform and the Utica Shale. Among these three units, the carbonate platform was assigned the highest hydraulic conductivity while the Utica Shale was assigned the lowest, creating a hydraulically isolated permeable unit in which water could travel faster than in the overlying and underlying geologic units. Slower flow above the carbonate platform results in an inverted distribution of groundwater residence times where older groundwater in the shale can be found above younger water in the carbonates (Figure 9).

Parametric study on faults

The Utica Shale formation (Figure 4) was recently targeted by the oil and gas industry as an unconventional natural gas reservoir. Public debate surrounding the environmental risks associated with such exploitation has raised many questions with respect to possible aquifer contamination (CCA 2014; Lavoie et al. 2014; Rivard et al. 2014; Lefebvre 2017).

While surface releases of fluids related to unconventional oil and gas operations (Vengosh et al. 2014) and well integrity (Nowamooz et al. 2015; Roy et al. 2016) pose the greatest threat to aquifers, geochemical evidence also seems to point to instances in which deep-seated fluids have migrated to the surface along natural preferential flow pathways, such as faults (Révész et al. 2010; Warner et al. 2012). Yet much debate still exists surrounding the possibility of shallow aquifer contamination through such pathways (Lefebvre 2017). While the question is often viewed from the perspective of the impact of fracking on fluid migration (Gassiat et al. 2013; Kissinger et al. 2013; Lange et al. 2013; Flewelling and Sharma 2014; Birdsell et al. 2015), natural mechanisms, such as regional flow which can contribute to the upward migration of natural formation fluids prior to exploitation, are less often explored.

Based on the observation of saline waters in valleys, Warner et al. (2012) hypothesised that the combination of high hydrodynamic pressure at depth under discharge areas (from regional flow, for instance) and the presence of higher permeability fracture zones can induce steep hydraulic gradients which could lead to migration of fluids between deeper strata and the surface. Indeed, this phenomenon seems to be occurring in the study area. Based on methane concentration observations, Pinti et al. (2013) suggested that the Yamaska Fault, an extension of the Jacques-Cartier River Fault (Figure 2), could act as a connection between deep formations and surface aquifers. Moreover, focusing their research on the Lotbinière study area, located at the northern end of the modelled cross-section (Figure 1), Lavoie et al. (2016) and Bordeleau et al. (2015, 2018) also found geochemical and physical evidence of upward flow near the Jacques Cartier River Fault, suggesting that it could act as a preferential pathway for migration of deep groundwater.

In this section, the potential role of a fault zone in a regional groundwater discharge area will be investigated through further numerical simulations. Since data gaps exist in the physical characterisation of the deep subsurface in the study area, this analysis does not attempt to deliver a definitive interpretation of the flow dynamics occurring near the Jacques-Cartier River Fault. The intent is rather to explore what impact different fault configurations would have on the hydraulic gradients and the distribution of mean groundwater ages near the fault zone as a basis for further studies.

Modelled fault configurations

In the context of this study, the fault zone is considered to be the volume of rock where the permeability has been altered as a result of fault-related deformations, while the protolith is considered the undeformed rock surrounding the fault zone (Bense et al. 2013). Resulting from tectonic activity, fault zones are complex hydraulic structures with various degrees of deformation. The fault core (FC) generally refers to the central zone of the fault which has undergone the most intense strain and most of the displacement, and as a result can have a lower permeability than the protolith (Bense et al. 2013). On the other hand, the damage zone (DZ), surrounding the fault core, is defined as the region which absorbed the remaining strain and as a result is marked by fracturing and an enhanced permeability (Bense et al. 2013).

The degree and nature of the rock deformation attributed to faulting depends on many intrinsic controls such as lithology, fault displacement and geometry, and deformation conditions (Caine et al. 1996) which are difficult to characterise. Furthermore, the definition of the fault’s permeability structure has also been shown to depend on the observation methods and geoscience discipline (Scibek et al. 2016).

Large-scale fault structures are often segmented and complex, exhibiting multiple fault cores and overlapping damage zones (Bense et al. 2013). In the current study, for the sake of simplicity and based on the conceptual model of Ladevèze (2017) proposed for the study area, the permeability structure of the Jacques-Cartier River fault zone is represented by a simplified conduit–barrier–conduit model (Figure 11). This model assumes the fault core is less permeable and the damage zone more permeable than the protolith.

Figure 11. Conceptual model of the Jacques-Cartier River fault zone. Vertical hydraulic conductivities shown correspond to scenario P of the parametric study on faults.

A number of fault characteristics such as the dip, the thickness and permeabilities of the fault zone and those of its core and its damage zone can play a role in controlling fluid flow around a fault (Bense et al. 2013). This study will focus on the impact of the permeability contrasts between the different parts of the fault zone by simulating a range of fault permeability structures (Table 4). Since the scale and nature of the hydrogeological numerical model herein is not suited for representing discrete fractures, the damage zone and fault core are represented using equivalent porous media. The geometry and dip of the fault are based on the geologic cross-section of Séjourné et al. (2013) (Figures 2 and 11) extending from the bedrock surface to the top of the crystalline formation of the Grenville Province. The thicknesses of the modelled fault core and the two damage zones each correspond to the width of one numerical grid element (50 m). For the different scenarios, the horizontal and vertical hydraulic conductivities of the fault core and damage zones are increased or decreased by up to three orders of magnitude with respect to the laterally adjacent elements corresponding to the protolith (Table 4), while retaining their original anisotropy. As the hydraulic conductivity of the Lorraine Group decreases with depth, so does that of the fault zone where it crosses the Lorraine Group. For comparison, tested scenarios also include conduit-only and barrier-only fault types.

Table 4. Simulated fault permeability scenarios.

Scenario	Fault core (FC) hydraulic conductivity contrast with protolith (P) in x and z directions	Damage zone (DZ) hydraulic conductivity contrast with protolith (P) in x and z directions	Fault zone permeability conceptual model
A	KFC = KP	KDZ = KP	Base case
B	KFC = KP x 0.1	KDZ = KP	Barrier
C	KFC = KP x 0.01	KDZ = KP	Barrier
D	KFC = KP x 0.001	KDZ = KP	Barrier
E	KFC = KP	KDZ = KP x 10	Conduit
F	KFC = KP x 0.1	KDZ = KP x 10	Conduit-barrier-conduit
G	KFC = KP x 0.01	KDZ = KP x 10	Conduit-barrier-conduit
H	KFC = KP x 0.001	KDZ = KP x 10	Conduit-barrier-conduit
I	KFC = KP	KDZ = KP x 100	Conduit
J	KFC = KP x 0.1	KDZ = KP x 100	Conduit-barrier-conduit
K	KFC = KP x 0.01	KDZ = KP x 100	Conduit-barrier-conduit
L	KFC = KP x 0.001	KDZ = KP x 100	Conduit-barrier-conduit
M	KFC = KP	KDZ = KP x 1000	Conduit
N	KFC = KP x 0.1	KDZ = KP x 1000	Conduit-barrier-conduit
O	KFC = KP x 0.01	KDZ = KP x 1000	Conduit-barrier-conduit
P	KFC = KP x 0.001	KDZ = KP x 1000	Conduit-barrier-conduit

Fault model results: Groundwater flow and age

Flow and age transport results for the various fault scenarios (Figures 12 and 13; see also figures C1 to C4 in annex C of Janos 2017) show that in the simulated flow system, a hydraulic conductivity contrast of at least two orders of magnitude between the protolith and the fault zone is required to induce significant diversion of the regional flow system. If the fault is acting as a barrier only, effects of the fault on the flow system begin to be observable when the fault zone hydraulic conductivity is two orders of magnitude less than that of the protolith but become truly significant with a difference of three orders of magnitude. Conversely, when the fault is acting as a conduit only, flow diversion is observable with a difference of only one order of magnitude in hydraulic conductivities between the fault and the protolith and, with a difference of two orders of magnitude, the conduit behaviour of the fault becomes clearly evident. This indicates that in the simulated system, the higher K of the damage zone exerts a greater influence on the flow system than the low K fault core (Figure 12).

Figure 12. Steady-state flow model showing hydraulic conductivities in the horizontal direction and streamlines for fault scenarios described in Table 4; vertical exaggeration is 1.25×.

Figure 13. Steady-state flow and age transport model showing mean groundwater age and streamlines for fault scenarios described in Table 4, vertical exaggeration is 1.25×.

From an age transport perspective (Figure 13), the presence of an effective barrier – that is, a fault core with a hydraulic conductivity of at least three orders of magnitude lower than that of the protolith – has the effect of lowering the active flow interface and decreasing the mean groundwater ages downstream of the fault (Figure 13d). By blocking regional flow, the barrier fault isolates the downstream flow system and local flow can penetrate deeper. Interestingly, an effective conduit-type fault, with a hydraulic conductivity of at least two orders of magnitude higher than the protolith, has a dual effect on the distribution of mean groundwater ages. For the regional flow system, the conduit-type fault zone is an effective preferential pathway, through which old water originating from long regional flow paths is concentrated and which may reach the surface (Figure 13i and m). However, in terms of the local flow system, the conduit-type fault also behaves somewhat as a drain, drawing water into the fault zone and thereby decreasing the mean groundwater ages and lowering the active flow interface upstream of the fault. This has the effect of displacing the old water plume from upstream to downstream of the fault (Figure 13i and m).

Combining a barrier-type fault core flanked with a conduit-type damage zone emphasises the effects of both the barrier and conduit zones (Figure 13p). On the upstream side of the fault zone, the damage zone channels regional flow, which is then effectively blocked by the fault core, allowing the plume of older water to reach closer to the surface before mixing with the local system. On the downstream side of the fault, the damage zone acts as a channel for local flow, allowing recharge to penetrate deeper into the subsurface and creating a larger mixing zone.

The most important outcome from this parametric study is that in the regional discharge area, the presence of a fault with a damage zone and/or a fault core with a significantly contrasting hydraulic conductivity with the protolith has the potential to divert and channel regional flow similar to a preferential pathway.

Nonetheless, this study represents a simplified model of the Jacques-Cartier fault and should only be viewed as a basis for further fault zone characterisation. To more accurately assess the impact of the Jacques-Cartier fault on the emergence of regional flow and its potential to act a pathway for contaminants related to fracking activities, further modelling should be performed taking into account the full complexity of fault zone hydrogeology. Such a model should refine the definition of the permeability structure of the fault zone and consider additional hydrodynamic processes which are relevant to deep basin flow such as density-dependent flow, pressure gradients in deep formations and multiphase flow.

Discussion

The aim of this study was to reinforce understanding of regional flow dynamics in the CA region through the analysis of characteristic time scales using geochemistry and numerical modelling. The conceptual model is based on Tóth’s (1999) topography-driven flow model and outlines the presence of local, intermediate and regional-scale flow systems.

The main regional aquifer was reported to be located in the upper fractured bedrock (Brun Koné 2013; Benoît et al. 2014; Lefebvre et al. 2015) and the numerical modelling presented herein confirms that the logarithmically decreasing hydraulic conductivity of the bedrock maintains active flow of groundwater in the top 40 m of the fractured rock. While simulated mean groundwater ages in the active flow zones are generally less than 100 years, travel times through the deeper local flow systems can reach more than 1 million years. The slightly shallower active and local flow systems in the St. Lawrence Lowlands, with respect to the Appalachian Highlands, highlight the effect of topographic relief on the local flow systems.

Rapid increases of the simulated mean groundwater ages with depth in the local flow systems appear to match the conceptual geochemical distribution model of the different water types (Figure 2). Water drawn from residential wells, which extend on average to depths of 60 m, would include young (G1-type) water from the active flow zone where flow rates are fastest, but would in some cases also include water from below the active flow zone – yielding the maximum radiocarbon age of 3500 years BP observed in this water group (Janos 2017). Furthermore, local flow systems in the study area can extend to maximum depths of approximately 200 m, and are likely affected by upward diffusion of old water from the regional flow system. A water composition gradient can therefore be expected within the same local flow system.

Moreover, the depth of 60 m in the model also roughly matches the zone where mean groundwater ages correspond to the retreat of the Champlain Sea. Water drawn from residential wells below the 10,000-year threshold will most likely be a mixture of salty water from the invasion episode and water that has recharged within the local flow system before the invasion. Given the relatively shallow local flow systems, in many cases it would also be highly probable for a residential well to intercept an intermediate or even the regional flow system.

Residential wells, however, are often installed as open boreholes in the bedrock aquifer, and thus the water they yield is a mixture of the different geochemical layers they cross. As shown by the flow model, groundwater flow velocities decrease drastically with depth; thus, contributions from the deeper layers of local flow systems and those of the regional system are expected to be less significant than those of the active flow layer.

The two water samples dated at 1.5 Ma by Bordeleau et al. (2015, 2018) were collected from a relatively shallow depth of 48 m near regional fault lines, suggesting regional faults such as the Jacques-Cartier River fault could act as preferential discharge pathways for regional flow. The parametric modelling study on fault zone permeability distributions showed that a difference of several orders of magnitude in hydraulic conductivities between the protolith, the damage zone and the fault core is required for the conduit/barrier effect of the fault to be significant. However, modelling also demonstrated that fault zones can indeed concentrate the input of regional flow to the active flow zone. Regional flow in the CA region is almost insignificant with respect to the active flow system and is therefore difficult to identify in the field. However, large-scale heterogeneities such as faults, which have the effect of concentrating flow, can enhance the detectability of regional flow in near-surface groundwater samples.

Conclusions

As part of the PACES aquifer characterisation programme, the conceptual and numerical flow and age transport models presented herein aimed to reinforce understanding of the regional flow system in the Chaudière-Appalaches region.

Based on a review of the regional geochemical context and the interpretation of radiocarbon ages of groundwater samples collected in the study area, groundwater flow was conceptualised as a nested flow system dominated by a series of shallow local flow systems. The 2D numerical groundwater flow model developed along a theoretical flow path confirms the presence of nested local, intermediate and regional flow systems.

Simulations of groundwater flow and age transport add to the evidence that groundwater flow in the CA region is mainly focused in the uppermost portion of the bedrock (< 40 m) where the fracture density is the highest, and suggest that near-surface fractures play an important role in controlling the flow. Below the active flow zone, mean groundwater ages in the local flow systems increase exponentially with depth. The effect of topography on the flow system is dominant in the Appalachian Highlands, which has well-developed intermediate flow systems. On the other hand, intermediate flow systems in the St. Lawrence Lowlands are almost non-existent and the interface between local and regional flow is relatively shallow (40–60 m).

While the fractured bedrock in the study area is well characterised near the surface (Crow and Ladevèze 2015; Lefebvre et al. 2015; Ladevèze et al. 2016), there is also evidence that major faults in the study area could likewise be associated with areas at a greater depth where the level of fracturing would be higher (Pinti et al. 2013; Lavoie et al. 2016; Bordeleau et al. 2018). Parametric modelling of a series of fault zone permeability distributions showed that a conduit, barrier or conduit–barrier–conduit-type fault could indeed act as a preferential pathway to regional groundwater flow and could explain the occurrence of geochemically evolved and old groundwater near such structures.

While the flow and mean age transport numerical model was useful for interpreting the general behaviour of the flow system, important simplifications were made due to data gaps. Additional field work focusing on better characterising the permeability structure near the Jacques-Cartier River Fault and the vertical distribution of groundwater salinity and mean age could help validate the presented conceptual model and provide further insight into the role of fault zones located in regional discharge areas. Further modelling, incorporating a variable density distribution and a transient flow system, would also deepen the understanding of the Chaudière-Appalaches regional flow dynamics.

Funding

The PACES-CA research project was funded primarily by the Quebec Environment Ministry (Ministère du Développement durable, de l’Environnement et de la Lutte contre les changements climatiques, MDDELCC), with contributions from many regional partners.

Ackowledgements

We acknowledge research support from the Natural Sciences and Engineering Research Council of Canada (NSERC), and from a Canada Research Chair in Quantitative Hydrogeology of Fractured Porous Media held by the second author. We are grateful for the assistance from the INRS PACES-CA team, notably Harold Vigneault, Jean-Marc Ballard, Marc-André Carrier, Marc Laurencelle, Xavier Malet, Annie Therrien, Châtelaine Beaudry and Guillaume Légaré Couture, for field work, technical assistance and advice. Christine Rivard, Geneviève Bordeleau and Pierre Ladevèze of the Geological Survey of Canada also provided generous support. We also thank P. Therrien of Université Laval who provided valuable computer technical support.