Stratigraphic sequence map for groundwater assessment and protection of unconsolidated aquifers: A case example in the Abitibi-Témiscamingue region, Québec, Canada

Abstract

Quantifying water resources at various scales is crucial for ensuring safe access to potable water for present and future generations. Fitting into this framework, this study presents a GIS-based approach aimed at allowing the evaluation of available groundwater resources within unconsolidated aquifers set in vast and heterogeneous shield regions. The approach was developed in a 19,397 km² region located in Abitibi-Témiscamingue (Québec, Canada) where unconsolidated aquifers are set in an irregular geologic framework owing to the rugged Canadian Shield topography and to the diversity of glacial and post-glacial events that shaped the landscape. Sparse and unevenly distributed stratigraphic borehole data were used for constructing a GIS-based model representing the total overburden thickness (unconsolidated sediments covering the bedrock), as well as the thickness and extent of the fine-grained deep-water sediments (silt-clay) of the Barlow-Ojibway proglacial Lake on a 100 m × 100 m mesh at the regional scale. The study suggests a new set of functions to map drift thickness in areas of highly discontinuous drift cover, using GIS tools. These data were used jointly with surficial deposits maps in order to develop a regional-scale two-dimensional model where 15 distinct stratigraphic sequences allow representation of the architecture of unconsolidated geological units. The approach allows simplifying data representation and ensuring the consistency between various regional-scale hydrogeological maps. The data is used to produce a stratigraphic sequence map for representing the architecture of aquifer– aquitard systems at the regional scale in a manner that is intelligible to non-specialists. The map and related data are first discussed for documenting the extent and volume of regional aquifers. The stratigraphic sequence map is subsequently discussed as a tool for supporting political decision makers dealing with issues related to groundwater resource assessment, protection and sustainable development.

La quantification des ressources hydriques est cruciale pour assurer un accès à l’eau potable pour les générations actuelles et futures. S’insérant dans ce contexte, cette étude présente une approche visant l’évaluation des ressources en eaux souterraines au sein d’aquifères non consolidés. L’approche a été développée pour une région de 19 397 km² située en Abitibi-Témiscamingue (Québec, Canada), dans un secteur où les aquifères non consolidés présentent une architecture irrégulière en raison de la topographie accidentée du Bouclier Canadien et de la diversité des événements glaciaires et postglaciaires ayant façonné le paysage. Des données de forages éparses et distribuées de façon hétérogène ont été utilisées pour la construction d’un modèle représentant l’épaisseur des dépôts meubles et des sédiments d’eau profonde (silt-argile) du lac proglaciaire Barlow-Ojibway sur un maillage de 100 m × 100 m. L’étude suggère de nouvelles fonctions pour cartographier à l’aide des SIG l’épaisseur de dépôts meubles dans des secteurs où la couverture est discontinue. Ces données ont été utilisées conjointement avec la cartographie des dépôts de surface afin de développer un modèle régional permettant la représentation à deux dimensions de l’architecture des unités géologiques non consolidées grâce à 15 séquences stratigraphiques. L’approche permet de simplifier la représentation des données et d’assurer la cohérence entre les différentes cartes hydrogéologiques régionales. Les données sont utilisées pour produire une carte des séquences stratigraphiques représentant l’architecture des systèmes aquifères–aquitards d’une manière intelligible pour les non-spécialistes. La carte et les données associées sont d’abord utilisées pour évaluer les limites et le volume des aquifères régionaux et ensuite présentées comme des outils permettant de soutenir les décideurs politiques qui traitent de questions liées à l’évaluation, à la protection et au développement durable des ressources en eau souterraine.

Introduction

Aquifers store the potable water resources for nearly half of the world’s population (World Water Assessment Program 2009). At the global scale, water resources (including groundwater) are increasingly stressed owing to population growth, human pressures and climate instabilities, among others (e.g. see Oki and Kanae 2006; Hoekstra and Mekonnen 2012). Under such circumstances, the protection and sustainable exploitation of groundwater resources at various scales stands as a critical issue in many countries. However, the development of land management strategies and political frameworks aimed at protecting groundwater must rely on robust hydrogeological data and expertise. This implies that groundwater protection inevitably depends on appropriate groundwater knowledge, as well as mutual knowledge transfer between the hydrogeological scientific community and political decision makers (e.g. see Alley and Leake 2004). There is a growing literature in hydrogeological sciences (e.g. see Niu et al. 2014), and scientists tend to focus on the use of tools and approaches that are increasing in complexity and sophistication. This represents a major advantage because it allows addressing hydrogeological problems at various scales with unprecedented efficiency and precision. However, knowledge transfer from scientists to political decision makers who are not experts in hydrogeology can become progressively more difficult if complex hydrogeological datasets need to be presented and explained. The key to this paradox may partly lie in the development of approaches for the simplified and robust representation of complex hydrogeological data.

The development of geological and hydrogeological models can be used to exemplify the issue outlined above. Models, both conceptual and computational, are central for understanding hydrogeological processes at various scales because they provide information about the extent, architecture and characteristics of geological units forming aquifers and aquitards and allow representation of flow processes. Conceptual models are useful because they allow simplified yet realistic representations of complex natural processes and provide a framework for sharing ideas with respect to the hydrogeological environment (Seibert and McDonnel 2002). However, conceptual models also suffer limitations when applied at the regional scale as they do not always incorporate the complexity of spatial variations in the geological framework. Computational models represent a valuable complement to conceptual models, as they allow representation of the hydrogeological environment within one-, two- and three-dimensional frameworks that reflect the interpretation of quantitative datasets at various scales. As of today, regional-scale computational geological models often rely on the combined use of GIS-based calculations and 3D modeling (Sharpe et al. 2007; Tremblay et al. 2010; Velasco et al. 2013; Van Hoesen 2014). Numerous challenges associated with the development of geological models have been addressed by the scientific community, and efficient approaches for interpolating data (Gao et al. 2006; Paulen et al. 2006; Boisvert and Parent 2008; Tremblay et al. 2010; Slaterry et al. 2011; MacCormack et al. 2013, Van Hoesen 2014) and developing hydrogeological interpretations (Ross et al. 2005; Jusseret et al. 2010; Lau et al. 2016) are increasingly documented. Various software programs (such as gOcad, GSI3D, Rockworks and ArcGIS) are used for constructing the models, and hydrogeological processes can be represented with unprecedented efficiency. Nevertheless, evaluating the architecture and extent of unconsolidated aquifers remains challenging, especially in areas where the bedrock topography is irregular and where quantitative data are sparse and unevenly distributed, and present variable quality/precision levels. In addition, illustrating the architecture of shallow unconsolidated aquifers over vast regions (i.e. > 10,000 km²) can be challenging owing to the contrast between the horizontal (tens of kilometers) and vertical (tens of meters) extent of models. Finally, there is a critical need to develop complementary approaches for allowing an intelligible representation of models in an effort to support decision makers with adequate hydrogeological knowledge. Such approaches could prove useful in improving exchanges and communications between scientists and political decision makers.

Fitting into this framework, the general objective of this study is to present an approach allowing the intelligible representation of 3D hydrogeological data as a tool for representing the architecture and extent of unconsolidated aquifers, discussing hydrogeological processes and supporting political decision makers dealing with issues related to groundwater resources assessment, protection and sustainable development. The specific objectives include (1) the development of GIS-based approaches for modeling the architecture of unconsolidated deposits using extensive borehole datasets in a region characterized by a rugged bedrock surface, and (2) the representation of this architecture on a 2D stratigraphic sequence map that is intelligible to non-specialists. The approaches are developed, applied and tested in a 19,397 km² area of the Canadian Shield in the Abitibi-Témiscamingue region (western Québec, Canada). The results are discussed to illustrate various applications related to the use of the aforementioned stratigraphic sequence map, to assess the state of the groundwater resource at the regional scale, to discuss regional-scale hydrogeological processes and to support decision makers in the development of land management strategies aimed at protecting aquifers. The study complements previous work conducted in other regions that focused on the joint evaluation of surface geology and overburden thickness (e.g. Soller 1993; Soller and Reheis 2004), and unconsolidated sediment stratigraphy and characteristics (e.g. Monahan et al. 2000). It provides new insights in relation to the quantitative assessment of aquifer architecture and extent in the context of the rugged Canadian Shield. The approach presented herein is best suited for regional-scale (i.e. 1000–10,000 km²) studies and is not intended to embrace the complexity of local-scale heterogeneities within the hydrogeological environment.

Study area

Covering 19,397 km², the study area (Figure 1) corresponds to the municipalities of the Abitibi-Témiscamingue region (western Québec, Canada) and is occupied by approximately 148,000 inhabitants (Ministère des Affaires municipales et de l'Occupation du territoire [MAMOT] 2016). It is estimated that more than 70% of this population relies on groundwater as its main source of drinking water, making it a critical resource at the regional scale (Cloutier et al. 2016). The municipal water systems supply about 76% of the population while the majority of unserved citizens use private groundwater wells. The main human activities affecting the landscape in the area include mining, sand and gravel extraction (Nadeau et al. 2015), agriculture and forestry, while many environmental liabilities related to mines and domestic waste are observed at the regional scale (Cloutier et al. 2016).

Figure 1. Study area. (a) Bedrock geology based on a compilation of geological maps at the 1:20,000 scale from the government of Québec, Système d'information géominière [SIGEOM] from the Ministère de l'Énergie et des Ressources naturelles du Québec [MERN]) (b) Surficial geology adapted from maps produced by the Geological Survey of Canada (Veillette 2004; Thibaudeau and Veillette, 2005; Paradis 2005, 2007; Veillette, 1986a, 1986b, 1987a, 1987b, 1987c); (c) elevation and continental water divide; (d) stratigraphic boreholes from geological and hydrogeological government databases.

Bedrock geology

The bedrock of the study region, based on the compilation of geological maps (at a 1: 20,000 scale) from the government of Québec (Système d'information géominière [SIGEOM] from the Ministère de l'Énergie et des Ressources naturelles du Québec [MERN]) (Figure 1(a)) mainly consists of Proterozoic and Archean intrusive, plutonic, volcanic and sedimentary rocks of the Superior Province (SP). A minor portion (307 km²) of the southern limit of the study area also overlaps with the Grenville Province. The region is intersected on an east–west axis by the Cadillac fault, a geological structure separating the Abitibi (ASP, north) and Pontiac (PSP, south) sub-provinces, both included in the SP. At the regional scale, the bedrock outcrops over approximately 10% of the territory and presents a rugged topography partly inherited from structural discontinuities (such as some regional faults) and differential erosion processes.

The ASP consists of a Neoarchean granite-greenstone geological domain hosting the Abitibi Greenstone Belt, one of the widest uninterrupted Neoarchean volcano-sedimentary belts of the planet (Hocq and Verpaelst 1994; Card and Poulsen 1998) and one of the most important mining regions of the world. The volcanic and sedimentary rocks of the ASP generally consist of strips elongated in the NW/SE orientation. The latter, characterized by a low-grade metamorphism (Card and Poulsen 1998), are intersected by numerous dykes and intruded by heterogeneously distributed felsic to mafic Archean plutonic rocks. The Cadillac fault defines the southern limit of the ASP within the study area.

South of the Cadillac fault, the PSP consists of a Neoarchean metasedimentary gneiss and granitoid geologic domain (Card and Poulsen 1998). Proterozoic sedimentary rocks of the Cobalt Group are identified mainly in the western portion of the study area, whereas Archean volcanic rocks are heterogeneously distributed.

In areas that are not served by municipal water systems, most people of the region withdraw groundwater from private wells installed in the fractured rock aquifers. The hydraulic properties of these fractured bedrock aquifers are highly variable at the regional scale owing to the diversity of host rocks and structural discontinuities found therein. Nevertheless, based on previous literature, Rouleau et al. (1999) suggested that water inflows in mines of the Canadian Shield can be abundant down to approximately 300 m below the bedrock surface owing to minor (e.g. joints) and major (e.g. faults and shear zones) structures. These authors also pointed out that sub-horizontal fractures are generally abundant within the first 75 m below the bedrock surface. Most of the private wells are set within this uppermost active layer.

Surface deposits

In North America, the Laurentide Ice Sheet reached its maximum extent approximately 20,000 years ago, during the last glaciation of the Wisconsin Episode (Quaternary Period). This last glaciation/deglaciation cycle represents the most recent geologic event that significantly contributed to the shaping of the landscape. In the southern portion of the study region, the ice began to retreat approximately 10,000 years ago (Veillette 1996) and the unconsolidated sediments are inherited from glacial and post-glacial events. Understanding the spatial distribution and architecture of these sediments is central for documenting the regional hydrogeological framework.

An exhaustive description of the regional surface geology can be found in the maps produced by the Geological Survey of Canada (GSC) (Veillette 1986a, 1986b, 1987a, 1987b, 1987c; Veillette 2004; Thibaudeau and Veillette 2005; Paradis 2005, 2007; hereafter referred to as ‘GSC maps’), as summarized in Figure 1(b). The synthetic stratigraphic column presented in Figure 2 (adapted from Veillette 1996) shows the main units of the region. For simplification purposes, these layers have been grouped into six distinct units that can present significant thickness contrasts according to local conditions.

Figure 2. Regional stratigraphic column (unconsolidated sediments) and associated simplifications used in the stratigraphic sequence map. For simplification purposes, the till unit is considered here to be composed of silt, sand and gravel.

Glacial deposits (till) represent the oldest unconsolidated unit of the region. This till unit, which is mainly inherited from crystalline rocks, consists of heterogeneous sediments composed of boulders, sand (≈ 60–70%), silt (≈ 20–25%) and clay (generally ≤ 5%) (Thibaudeau and Veillette 2005). At the regional scale, the till unconformably overlies the Precambrian bedrock. It outcrops over approximately 14% of the territory, mainly in the vicinity of bedrock outcrops and in sectors of higher elevation where younger deposits are absent. Buried till patches are also likely present in different bedrock depressions, beneath younger deposits. Nevertheless, the identification of buried till deposits from pre-existing borehole data is often challenging owing to its discontinuous distribution and its similarity (in terms of matrix grain size distribution) to some of the overlying sediments. The extent of the till unit beneath younger sediments therefore remains poorly documented within the study region.

Glaciofluvial sediments were deposited by glacier meltwaters. In the region, they form eskers that are characterized by a coarse granular core flanked by subaqueous fans developed at the emergence of subglacial meltwater at the ice front, within proglacial lake Barlow‐Ojibway. These eskers form a vast network converging toward a large interlobate landform (Harricana Moraine, Figure 1(b)) oriented roughly north–south. The Harricana-Lake McConnell glaciofluvial system was the expression used by Veillette (1986c) to refer to the southern component of this extensive and complex landform (Nadeau et al. 2015). The moraine shares similarities with most of the eskers of the regions in terms of depositional processes and grain size distribution (e.g. see Veillette 1986c, 1988; Nadeau et al. 2015). The latter generally consist of coarse grained complexes (sand and/or gravel matrix according to the gradistat scale (Blott and Pye 2001)) altogether outcropping over 4.6% of the territory and often defining elongated sinuous high-relief features within the landscape owing to the preferential accumulation of coarse-grained sediments. Some of these formations are recognized as the most productive aquifers of the study region (Cloutier et al. 2016). Subaqueous outwash sediments mainly consisting of sand are also widespread between the eskers and moraines of the region but are most often covered by younger sediments.

Fine-grained deep-water glaciolacustrine sediments (silt–clay) of the Barlow-Ojibway Lake compose the featureless clay plain occupying 41.5% of the study area. The hydrogeological characteristics of this unit remain poorly documented at the regional scale, mainly because it is rarely exploited for groundwater supply owing to its grain size characteristics and to its generally assumed low hydraulic conductivity (e.g. see Cloutier et al. 2013, 2015, 2016).

Sublittoral and deltaic sediments are associated with the progressive retreat of glaciolacustrine waters. Sublittoral sediments mainly consist of broad surfaces of sand and gravelly sand overlying fine-grained deep-water glaciolacustrine sediments (Nadeau et al. 2015). Deltas are mainly associated with eskers and occupy a transition zone between glaciofluvial and glaciolacustrine deposits. Altogether, these sediments cover approximately 10.7% of the study area.

Alluvium deposits represent the uppermost inorganic unit of the regional stratigraphic column. These deposits mainly consist of gravelly sand, sand, silty sand and clayey silt that were deposited in flood plains. This unit rarely exceeds a few meters in thickness and occupies less than 2% of the study area.

Organic deposits cover 16% of the study area, mainly on esker/moraine flanks and on the fine-grained glaciolacustrine sediments. Despite their widespread occurrence, organic deposits rarely exceed 5 m in thickness (e.g. see Cloutier et al. 2015).

Geomorphology, hydrology, climate and vegetation

Nadeau et al. (2015, 537) described the geomorphology of the region as ‘a featureless clay plain broken up locally by wave-washed bedrock knobs, higher hills covered by a thin cover of till and glaciofluvial deposits, with eskers and moraines rising above the surrounding terrain’. Surface elevations (Figure 1(c)) range between 174 and 567 meters above sea level (asl) and the highest points correspond to sparsely distributed bedrock outcrops, whereas several elongated and sinuous high-relief features correspond to the presence of eskers and moraines. Despite its relatively flat terrain, the study area encompasses the continental water divide between waters flowing into the James Bay (north) and the St. Lawrence River (south) basins (Figure 1(c)). Overall, 18,488 ponds (< 0.1 km²) and 489 lakes (> 0.1 km²) are identified within the area, and wetlands (swamps and peatlands) altogether cover more than 4818 km² at the regional scale (Cloutier et al. 2016). The vegetation of the area is dominated by coniferous species (spruce, fir and pine) in the north, whereas a transition to deciduous species (mainly birch) is observed southwards. At the Mont-Brun meteorological station (Figure 1(c); Environment Canada station #7085106) annual rainfall and snowfall average approximately 705 mm and 281 mm, respectively, whereas daily average air temperatures range between approximately −18 °C in January and 17 °C in July (Environment Canada, 1981 to 2010 Canadian Climate Normals station data).

Methods

The approaches described herein ultimately aim to create a stratigraphic sequence map consisting of a regional-scale 100 m × 100 m raster that allows an integrated representation of surface deposit and unconsolidated sediment architecture. Although a proper knowledge of unconsolidated sediment thicknesses is required by the approach, the stratigraphic sequence map consists of a qualitative representation of the geological environment as the x–y information from typical surface deposits maps is respected whereas the succession of stratigraphic layers along the z-axis is illustrated for each cell of the raster map, without representing their thicknesses. The concept supporting the approach is illustrated in Figure 3, where the main data types and methodological steps required to construct the stratigraphic sequence map are shown. The latter are discussed separately below.

Figure 3. Conceptual representation of the approach for constructing the stratigraphic sequence map.

Digital elevation model

A regional scale digital elevation model (DEM) was constructed on a 10 m × 10 m grid using data from the Banque de données topographiques du Québec (BDTQ) at a 1:20,000 scale and accounting for features of the structured hydrographic network (data from the Centre d’Expertise Hydrique du Québec). The DEM was interpolated using the Topo to raster tool in ArcGIS 10.2 which uses the ANUDEM elevation gridding procedure (Hutchinson et al. 2011) and allows representation of the hydrographic network. In order to reduce gaps at the junction of cartographic sheets, interpolations were conducted separately for groups of watersheds using a 500-m external buffer. The different groups of watersheds were subsequently regrouped and average values were used for overlapping areas corresponding to the 500-m external buffers. The precision of the DEM is assessed by comparing pixel elevation values to reference data for 272 selected sites of the provincial geodetic network that were selected in order to cover a wide range of elevations and to avoid anthropic structures such as bridges and buildings. The comparison revealed maximal error of 9.6 m and a regional mean absolute error of 2.4 m. A coarser (i.e. 100 m × 100 m mesh size) DEM was created based on average values of the 10 m × 10 m DEM for facilitating regional-scale calculations. Both DEMs are employed in the approaches described below.

Surface roughness measurements

Surface roughness calculations are needed to assess the spatial representativeness of borehole data and must be considered when interpreting regional-scale maps of unconsolidated sediment thickness. Surface roughness is described herein as the absolute amplitude of changes in elevation that can be measured within predefined circle diameters (30, 50, 90, 250 and 490 m) for a given group of outcropping geological units. The DEM with 10 m × 10 m cells was used for calculating surface roughness (SR) for selected groups of outcropping geological units. For simplification purposes, three groups were identified within the regional geological units:

•	Group 1: Bedrock outcrops and shallow till patches;
•	Group 2: Glaciofluvial, sublittoral and deltaic deposits;
•	Group 3: Fine-grained deep-water glaciolacustrine sediments.

The three groups of geological units define distinct trends of increasing SR vs. distance, as graphically shown in Figure 4. Overall, Group 1 presents the highest SR. For instance, Table 1 suggests that average variations in bedrock topography reaching 28.4 m are likely to be observed within a circle of 490 m in diameter.

Figure 4. Average range of surface elevation (referring to surface roughness, SR) for groups of geological units within predefined diameters.

Table 1. Average range (in meters) of surface elevation (surface roughness, SR) and associated standard deviations (σ) within a predefined circle diameter (distance) for groups of geological units.

Distance (m)	Group 1			Group 2			Group 3
	SR	σ	SR+2σ	SR	σ	SR+2σ	SR	σ	SR+2σ
30	1.4	1.4	4.26	0.6	0.8	2.1	0.55	0.71	1.97
50	2.8	2.8	8.32	1.2	1.4	4.1	1.07	1.33	3.73
90	5.7	5.2	16.11	2.4	2.7	7.7	2.11	2.37	6.85
250	16.1	11.6	39.34	6.7	5.9	18.6	5.41	4.41	14.23
490	28.4	17.2	62.73	12.9	9.3	31.4	9.12	5.67	20.46

Group 1: Bedrock outcrops and shallow till patches; Group 2: Glaciofluvial, sublittoral and deltaic deposits; Group 3: Fine-grained deep-water glaciolacustrine sediments.

Unconsolidated sediment thickness

Table 2 identifies the data sources and summarizes the GIS-based preliminary treatments that were conducted for preparing the dataset, whereas Figure 1(d) illustrates the spatial distribution of available data. It is worth mentioning that the diamond-drilling boreholes (n = 30,195) only provide total overburden thickness, without further stratigraphic information with respect to surface deposits. The wells from the Hydrogeological Information System (HIS, as reported in Table 2) only present simplified information with respect to stratigraphy. These data allow documenting the thickness of organic, fine-grained and coarse-grained sediments, without specific geologic information. To the authors’ knowledge, there is currently no quantitative evaluation of the error associated with the unconsolidated sediment thickness (UST) provided with the different borehole databases. The eight main steps related to the evaluation of UST are summarized in Figure 5(a) and described separately below. The procedure illustrated therein was developed in order to account for the main constraints related to the regional framework:

(1)	The available borehole data are unevenly distributed;
(2)	The Precambrian Shield is rugged and presents significant SR;
(3)	The available borehole data are unlikely to allow a robust representation of the elongated hollows and crests associated with regional faults and glaciofluvial formations.

Table 2. Data sources and number of boreholes (n) for evaluating unconsolidated sediment thickness (UST).

Data types	Data sources	n	Preliminary GIS data treatments
Bedrock outcrops	- Government database^1	132,651	A value of zero is attributed to overburden thickness for each point-specific bedrock outcrop. In the case of bedrock/thin till outcrops defining polygons on the GSC maps, points of zero overburden thickness are created at vertices of the contour and in the center of corresponding 100 m × 100 m cells.
	- Point-specific outcrops from GSC maps^2	27,912
	- Surface and contours of polygons corresponding to bedrock and shallow till outcrops in the GSC maps^2	717,345
Boreholes	- Wells from the Hydrogeological Information System (HIS) government database^3	3918	A systematic approach is applied in order to subtract inconsistent data: - Boreholes identified within bedrock/thin till outcrops of the GSC maps; - Boreholes revealing a zero unconsolidated sediment thickness; - Duplicate data; - Boreholes suggesting unrealistic unconsolidated thicknesses (i.e. >100 m); The position (x, y, z) of data are corrected based on borehole dip and direction.
	- Mining diamond drillings^1	30,195
	- Overburden drilling management (ODM) from government database^1	1282

¹ Ministère de l’Énergie et des Ressources naturelles du Québec.

² Geological Survey of Canada.

³ Ministère du Développement durable, de l’Environnement et de la Lutte contre les changements climatiques du Québec.

Figure 5. GIS-based procedures for evaluating (a) the unconsolidated sediment thickness (UST) and (b) the fine-grained glaciolacustrine unit thickness (FGGUT). DNBO: distance to nearest bedrock outcrop; MF: multiplication factor; MUST: maximum unconsolidated sediment thickness; DEM: digital elevation model; DNL: distance to nearest limit of the clay plain.

Calculation of regional trends

A straightforward interpolation between borehole and outcrop data would not allow realistic regional-scale representations of UST owing to the rugged Precambrian Shield topography. The main issue is related to the lack of borehole data for documenting UST between bedrock outcrops, causing underestimates in UST when straightforward interpolation is used. The key hypothesis supporting the alternative approach proposed herein is that on average, UST is likely to increase with increasing distance from a bedrock outcrop (as suggested in Boisvert and Parent 2008). This hypothesis can be tested by grouping boreholes in classes according to the omnidirectional distance to the nearest bedrock outcrop (DNBO) and calculating average UST for each class (Figure 5(a) – steps 1–2). The regional-scale trend that is obtained by applying this procedure is shown in Figure 6, where the data suggest that average UST increases asymptotically with increasing DNBO and that a natural logarithmic function can be used for representing the data for DNBO ≤ 1900 m:

Figure 6. Regional trend of average unconsolidated sediment thickness (UST, solid line) vs. distance to nearest bedrock outcrop (DNBO) and maximum plausible unconsolidated sediment thickness (MUST) equal to the average UST, plus surface roughness of bedrock (SR), plus twice the standard deviation on the SR ( in m), both with reference to the left-side y-axis. The number of available borehole data (dashed line) with respect to each class of DNBO is shown, with reference to the right-side y-axis.

(1)

where and are both expressed in meters. Beyond the 1900-m cut-off, the paucity of available data precludes the use of robust regional averages. Nevertheless, the data suggest that USTs do not systematically increase beyond that limit. Equation (1) is used to produce a preliminary/intermediate UST raster based on the calculated DNBO for each 100 m × 100 m cell of the regional grid (Figure 5(a) – step 3).

Calculation of local trends

Because the calculations presented above are based on regional-scale averages, they do not allow a robust representation of local heterogeneities in UST. As an example, boreholes located at a distance of 100 m from an outcrop have a calculated regional average UST of 9.84 m, although UST values ranging from 0 to 99 m are identified among these data. This suggests that an envelope (delineated by at least two functions in an UST vs. DNBO graph) rather than a single function should be used for estimating UST at the regional scale. Using a multiplication factor (MF) for adjusting the outputs of Equation (1) allows for better representation of local-scale UST from the borehole dataset:(2)

where is the UST associated with a specific borehole and is the UST deduced from Equation (1) for the same borehole. It is assumed here that UST values calculated from a direct multiplication of values (Equation 1) by the MF values (Equation 2) are likely to represent a significant overestimate of true UST in sectors where extreme (high) MF values are calculated and/or borehole data are inaccurate. A third function is therefore proposed for defining the maximum plausible UST vs. DNBO trend. The latter is established based on Equation (1) and on the SR defined for Group 1. It is proposed here that for a given DNBO class, the maximum plausible unconsolidated sediment thickness (MUST) is equal to the result associated with Equation (1), plus SR, plus twice the standard deviation on the SR ( in m) (also see Table 1). The MUST calculated using this approach increases asymptotically with increasing DNBO (Figure 6), and a natural logarithmic function can be used for representing the data:

(3)

The USTs from boreholes that do not allow respecting the MUST criteria (Equation 3) are removed from the dataset in order to reduce UST overestimates (Figure 5(a) – step 4). The remaining boreholes (n = 34,743) are subsequently interpolated by inverse distance weighting to produce a raster map (Figure 5(a) – step 5). The interpolated MF is used to multiply the UST values evaluated from Equation (1) to produce an estimate of UST () that accounts for local borehole data and/or bedrock SR (Equation 4; Figure 5(a) – step 6):

(4)

Based on Equation (4), the UST₂ raster represents a corrected version of the UST₁ raster. This procedure allows better representation of local trends in UST vs. DNBO values.

Correction for the geomorphology of glaciofluvial deposits

A regional-scale preliminary raster of bedrock topography can be constructed by subtracting the UST₂ values (Equation 4) from the values of the DEM. Nevertheless, this bedrock topography raster presents a major inadequacy, as crests in bedrock topography are systematically associated with the large glaciofluvial formations of the region. This reveals that an additional correction needs to be applied in order to account for the geomorphological features associated with glaciofluvial sediments (Figure 5(a) – step 7). Because most of these formations consist of elongated features characterized by preferential sediment accumulation where borehole data are often unavailable, they are not well-represented by the previously described calculations, which are more appropriate for representing the infilling of bedrock depressions by sediments. This shortcoming is addressed by projecting the vertices of the polygons delineating glaciofluvial formations on the preliminary bedrock topography raster in order to create a triangular irregular surface (TIN) that is used to truncate the previously described artifact and to add the corresponding value to the UST. This allows an improved representation of the preferential sediment accumulation associated with long and narrow glaciofluvial formations (UST_final in Figure 5(a) – step 8).

Fine-grained glaciolacustrine unit thickness

The data sources and preliminary GIS-based procedures that were realized for preparing the dataset in view of evaluating the extent of the fine-grained glaciolacustrine unit (FGGU) are summarized in Table 3. Following data preparation, a reference surface is created by interpolation for representing the maximum elevation reached by the FGGU at the regional scale (Figure 7):

(1)	Points are created at the centroid of segments defining the contact between the FGGU and underlying units (e.g. bedrock and till) in the GSC maps. The elevations of these points are extracted from the DEM and interpolated by inverse distance weighting to produce a surface (raster map) representing the maximum elevation of the clay layer at the regional scale (Figure 7 – step 1);
(2)	The aforementioned surface is suppressed in sectors where its elevation is lower than that of the DEM (Figure 7 – step 2). In some cases, especially in the vicinity of esker flanks, the elevation of FGGU will be higher than for younger deposits, suggesting that the FGGU is found under the more recent deposits;
(3)	The surface is expanded to account for sectors where the FGGU is identified at the surface in the GSC maps (Figure 7 – step 3) to define a polygon representing the clay plain limit.

Table 3. Data sources and number of boreholes (n) for evaluating the fine-grained glaciolacustrine unit thickness (FGGUT).

Data types	Data sources	n	Preliminary GIS data treatments
Stratigraphic drilling	- Wells from Hydrogeological Information System (HIS) government database^1	3061	Selection of boreholes for which (1) the FGGU is identified as the dominant facies for at least one layer and (2) the total unconsolidated sediment thickness is known. Calculation of the total FGGUT for each of the boreholes selected above. Suppression of duplicate data.
	- Overburden drilling management (ODM) from government database^2	989
	- Technical reports and work by Cloutier et al. (2013, 2015)	199

¹ Ministère de l’Énergie et des Ressources naturelles du Québec.

² Ministère du Développement durable, de l’Environnement et de la Lutte contre les changements climatiques du Québec.

Figure 7. Conceptual GIS-based method to define the reference surface delineating the maximum elevation of the fine-grained glaciolacustrine unit (FGGU). The surface is obtained by interpolating the highest elevations reached by the FGGU.

Calculation of regional trends

Following the construction of the polygon representing the extent of the FGGU, the available stratigraphic data with respect to the FGGU (Table 3) are grouped in classes according to the distance to the nearest limit (DNL; in a similar fashion to what was proposed for UST evaluation) to evaluate an average FGGU thicknesses (FGGUT) function (Figure 5(b) – steps 1–2). The data reported in Figure

Figure 8. Regional trend of average fine-grained glaciolacustrine unit thickness (FGGUT, solid line) vs. distance to the nearest limit of reference plan (DNL), with reference to the left-side y-axis. The number of available borehole data (dashed line) for each class of DNL is shown, with reference to the right-side y-axis.

8 suggest that average FGGUT increase with increasing DNL and that a linear function can be used for representing the data for distances ≤ 1300 m:

(5)

where and are both expressed in meters. Beyond the 1300-m cut-off, the paucity of available data precludes the use of robust regional averages. Equation (5) is used to prepare a preliminary raster of FGGUT (Figure 5(b) – step 3).

Calculation of local trends

A multiplication factor (MF_FGGUT) is used for adjusting the outputs of Equation (5) in order to allow better representation of local-scale FGGUT vs. DNL trends:

(6)

where is the FGGUT associated with a specific borehole, and is the FGGUT deduced from Equation (5) for the same borehole (Figure 5(b) – step 4). The MF_FGGUT (Figure 5(b) – step 5) are interpolated by inverse distance weighting and subsequently used to multiply the FGGUT values evaluated from Equation (5) to produce a regional-scale representation of FGGUT () that accounts for local features (Equation 7; Figure 5(b) – step 6):

(7)

A final correction is made to keep the minimum value between FGGUT₂ (Equation 7) and the UST_Final to obtain the final FGGUT layer, which cannot exceed the total UST (Figure 5(b) – step 7).

Creation of a stratigraphic sequence map

The stratigraphic sequence map consists of a rasterized surface deposit map where each cell is associated with its own stratigraphic column (as initially conceptualized in Figure 2). The estimated thickness of UST and FGGUT are used for constructing the regional framework with 15 different stratigraphic sequences. The logical approach, and the criteria used for building GIS-based queries and partitioning the environment according to these 15 sequences, are reported in Table 4. The aforementioned criteria are based on the following simplifications and assumptions:

(1)	The extent of the till layer is discontinuous at the scale of the region and this unit can likely be found within stratigraphic sequences 2 to 15, although the approach developed herein does not allow representing the thickness of this unit;
(2)	In areas where the FGGUT is overlain by sublittoral sand (mainly on esker flanks, stratigraphic sequences #6–7), the maximum thickness of the sand unit is restricted to 5 m, although this unit can reach more significant thicknesses within the region;
(3)	A thickness of 5 m is allowed for organic deposits (stratigraphic sequences #11–12), although peatlands of the region do not systematically reach this vertical extent;
(4)	When constructing the stratigraphic sequence map, the information associated with the GSC maps (surface geology) prevails over the other data sources and/or calculations. As an example, if the calculations suggest that the UST = 0 but the surface geology map allows identifying sediments at the surface, a thin layer of unconsolidated sediments is assumed to be present and is represented in the stratigraphic sequence map.

Table 4. Description of the criteria for defining the 15 regional stratigraphic sequences.

SS#	Criteria for defining the stratigraphic sequences
	Surface deposits (based on GSC maps)	UST (m)	FGGUT (m)	UST − FGGUT (m)	Legend scheme
1	Bedrock and shallow till	= 0	= 0	= 0
2	Till	> 0	= 0	> 0
3	Glaciofluvial	> 0	= 0	> 0
4	Glaciolacustrine deep-water	> 0	> 0	> 0
5	Glaciolacustrine deep-water	> 0	> 0	= 0
6	Glaciolacustrine littoral and deltaic	> 0	> 0	≥ 5
7	Glaciolacustrine littoral and deltaic	> 0	> 0	< 5
8	Glaciolacustrine littoral and deltaic	> 0	= 0	> 0
9	Organic	> 0	> 0	> 0
10	Organic	> 0	> 0	= 0
11	Organic	≥ 5	= 0	≥ 5
12	Organic	< 5	= 0	< 5
13	Alluvium	> 0	> 0	> 0
14	Alluvium	> 0	> 0	= 0
15	Alluvium	> 0	= 0	> 0

SS#: Stratigraphic sequence number; GSC: geological survey of Canada; UST: unconsolidated sediment thickness; FGGUT: fine-grained glaciolacustrine unit thickness; UST − FGGUT: difference between UST and FGGUT values.

Aquifer delineation

The fine-grained glaciolacustrine unit (hereafter referred to as the clay aquitard) is assumed to present a hydraulic conductivity that is significantly lower than that of the other inorganic unconsolidated units (see Cloutier et al. 2013, 2015, 2016). Based on this assumption, the sand and gravel deposits found beneath the clay aquitard (stratigaphic sequences #4, 6, 9 and 13) are systematically categorized as confined aquifers. Numerous field observations (Cloutier et al. 2013, 2015) and the widespread occurrence of point-specific and diffuse springs found at the margin of glaciofluvial formations, at the contact with the clay unit (Nadeau et al. 2015), support the hypothesis that the sediments found beneath the regional aquitard are most often water saturated. The volumes of confined aquifers can therefore be estimated solely on the basis of geological data, without having to provide an extensive interpretation of groundwater levels.

Although some sand and gravel sediments found above the clay aquitard are most likely to store groundwater (especially the uppermost inorganic units in stratigaphic sequences #6 and 7; Table 4), these formations are generally limited in terms of volume within the study region owing to their shallow thicknesses (generally ≤ 5 m) and are rarely exploited for groundwater supply.

The main unconfined granular aquifers are assumed to be restricted to sectors where the clay aquitard is absent and mainly associated with stratigaphic sequences #3, 8 and 11. Estimating groundwater levels within these formations is challenging, owing to the lack of available data at the regional scale. Nevertheless, based on the above-mentioned hypothesis stating that the sediments found beneath the regional aquitard are most often water saturated, it is postulated here that the interpolated maximum elevation of the clay layer (FGGU) allows creating a reference surface that provides a realistic estimate of the minimum plausible groundwater level within unconfined granular aquifers (as depicted in Figure 7). The aforementioned approach is consistent with previous work by Nadeau et al. (2015), who suggested that the elevation of groundwater within the eskers and moraines of the region is imposed by the elevation of the clay aquitard and/or bedrock outcrops delineating their extent and can be represented by an interpolated reference surface. The data from Nadeau et al. (2015) suggest that although this reference surface is approximate, it allows the representation of groundwater levels within eskers in a manner that is consistent with (1) the available data from groundwater wells, (2) the regional geological framework, and (3) the spatial distribution of groundwater springs and diffuse seepage zones identified within the area (Cloutier et al. 2013, 2015, 2016). In that sense, it is assumed that the surface representing the maximum elevation of the FGGU layer at the regional scale allows a realistic estimate of the position of the groundwater table. The thickness of unconfined granular aquifers can thus be estimated by subtracting the elevation of the bedrock surface from that of the reference surface.

Results

The results associated with the calculation of UST (Figure 9(a)) and of FGGUT (Figure 9(b)), when combined with the surface deposits map (Figure 2(b)), allow construction of the stratigraphic sequence map (Figure 9(c)) and realization of a quantitative assessment of unconsolidated aquifer extent at the regional scale (Figure 9(d)).

Figure 9. Main regional-scale outputs of the GIS-based calculations. (a) Unconsolidated sediment thickness (UST); (b) fine-grained glaciolacustrine unit thickness (FGGUT); (c) stratigraphic sequence map; (d) calculated aquifer thicknesses. The black rectangle corresponds to the area covered in Figure 10 which provides local-scale data for a sector located west of the city of Amos.

As a complement to the regional maps (Figure 9), Figure 10 shows the UST (Figure 10(a)), FGGUT (Figure 10(b)) and surface geology (Figure 10(c)) data at the local scale, with a cross section (A–A’) showing the corresponding morphology of the bedrock and overlying sediments (Figure 10(d)). The numbers shown in the cross section of Figure 10(d) correspond to the stratigraphic sequences, as identified in Table 4, consistent with the 2D plan view of the local stratigraphic sequence map (Figure 10(e)). Although the approach presented herein is mainly applicable at the regional scale, Figure 10 is used to exemplify the consistency between the UST, FGGUT and surface geology maps, and the corresponding stratigraphic sequence map.

Figure 10. Local-scale illustration of the consistency between the maps and datasets supporting the construction of the stratigraphic sequence map. (a) Unconsolidated sediment thickness (UST). (b) Fine-grained glaciolacustrine unit thickness (FGGUT). (c) Surface geology map. (d) A–A` cross section showing the bedrock (red) and overlying coarse-grained (orange) and fine-grained (purple) sediments. The numbers reported above the cross section correspond to the different stratigraphic sequences, as identified in Table 4. (e) Corresponding local-scale stratigraphic sequence map with numbers corresponding to the different stratigraphic sequences identified in Table 4.

The extent of unconsolidated units: Aquifer delineation and resource assessment

Unconsolidated sediment thicknesses: Regional trends

The results reported in Table 5 suggest that regional UST values can reach up to 91 m, for a regional average of 11.2 m and a total volume of 198 km³. The FGGUT values reach up to 66 m, for a regional average value of 9.3 m and a total regional volume of 93.1 km³ (accounting for 47% of the total unconsolidated sediments). The regional maps (Figure 9(a–c)) suggest that areas of preferential sedimentation can be divided into two groups:

(1)

The main eskers and moraines (Figure 1(b); stratigraphic sequence #3, Figure 9(c)) of the region define sinuous crests within the regional landscape and correspond to areas of positive sediment accumulation. The spatial distribution of eskers and moraines is most likely related to a combination of glacial processes that occurred within the Laurentide Ice Sheet and of sediment accumulation processes at the ice margin within the Barlow-Ojibway proglacial Lake. Regionally, eskers and moraines are oriented roughly along the north–south axis, and their thicknesses seem to be related to their widths.

(2)

The sectors associated with bedrock depressions correspond to areas where sediments have preferentially accumulated owing to limited erosion and/or gravitational processes. Narrow bedrock depressions/hollows found at various elevations can be associated with sectors where the till deposits have been preserved from post-glacial erosion. The vast sectors of low bedrock elevations associated with lakes Abitibi and Témiscamingue correspond to the areas characterized by the greatest UST values.

The greatest FGGUT values (Figure 9(b)) are located in sectors corresponding to the highest UST and associated with vast bedrock depressions. These sectors most likely corresponded to the deepest areas of the Barlow-Ojibway proglacial Lake, where a limited energy within the water column allowed the preferential accumulation of fine-grained sediments. These sectors also correspond to areas where the remnants of lakes Barlow (roughly located south of the Cadillac fault) and Ojibway (roughly located north of the Cadillac Fault) persisted over the longest period between ice retreat and lake drainage (Veillette 1996; Roy et al. 2015), allowing more time for sediment accumulation.

Table 5. Descriptive statistics and volume estimates in relation to the 15 regional stratigraphic sequences.

Stratigrpahic sequence	Area (km^2)	% of study area ^i	Elevation					UST (m)				FGGUT					Confined granular aquifers					Unconfined granular sediments				Unconfined granular aquifer
			Min (m)	Max (m)	Range (m)	Mean (m)	SD (m)	Max (m)	Mean (m)	SD (m)	Volume (km^3)		Max (m)	Mean (m)	SD (m)	Volume (km^3)	Max (m)	Mean (m)	SD (m)	Volume (km^3)	Groundwater volume ii (km^3)	Max (m)	Mean (m)	SD (m)	Sediment volume (km^3)	Max (m)	Mean (m)	SD (m)	Volume (km^3)	Groundwater volume ii (km^3)
0	1774
1	3101	17.6	180	567	387	315	34
2	1313	7.5	176	556	380	333	32	67	8	7	10											67	10	6	11
3	808	4.6	175	426	251	326	29	91	23	13	18											91	23	13	19	90	15	10	6	2
4	4844	27.5	180	381	201	296	22	87	17	7	80		66	10	5	49	71	7	6	32	11
5	2469	14.0	180	377	197	293	25	64	7	6	17		64	7	6	17
6	361	2.1	180	366	186	304	31	78	23	10	8		65	9	6	3	63	8	8	3	1
7	351	2.0	180	364	184	285	41	69	13	9	4		65	11	9	4
8	1184	6.7	175	428	253	328	22	74	12	9	14											74	13	8	14	66	11	8	6	2
9	1319	7.5	180	370	190	303	18	83	17	7	23		58	10	4	13	64	8	6	10	4
10	366	2.1	180	357	178	297	23	59	11	7	4		59	11	7	4
11	905	5.1	229	428	199	324	13	56	15	7	13											56	15	7	13	56	11	7	8	3
12	234	1.3	185	439	255	314	21	5	1	2	-											5	4	1	-	5	2	1	-	-
13	234	1.3	180	346	166	294	23	64	17	8	4		47	10	5	2	57	7	6	2	1
14	83	0.5	180	335	154	291	24	55	10	7	1		55	10	7	1
15	53	0.3	183	457	274	313	25	42	6	7	-											42	11	6	-	38	9	7	-	-
Total	19,397	100	175	567	392			91	11		198		66	9		93	71	7		46	16	91	15		57	90	7		21	8

ⁱ The areas corresponding to lakes are excluded.

ⁱⁱ Groundwater volumes are calculated assuming a porosity of 35%. UST: Unconsolidated sediment thickness. FGGUT: Fine-grained glaciolacustrine unit thickness. SD: Standard deviation.

Overall, at the regional scale, the UST values calculated for each 100 m × 100 m parcel of the region illustrate the link between bedrock surface elevation and sediment accumulation (Figure 11). Two trends are shown in Figure 11, where the average UST values gradually increase as the bedrock surface elevations gradually decrease, from the highest portions of the study area that are generally associated with bedrock outcrops (i.e. sites where UST = 0 m) toward the low-lying sectors associated with Lake Abitibi (northern portion of the study region; remnant of Lake Ojibway) and Lake Témiscamingue (southern portion of the study region; remnant of Lake Barlow). The data presented in Figure 11 also suggest that in sectors where the bedrock surface is at an elevation greater than 300 m, the thickest unconsolidated sediments correspond to glaciofluvial deposits. Another key observation is that the difference between total average UST and glaciofluvial UST is generally smallest in the vicinity of Lake Abitibi. This suggests that the eskers found in the northwestern portion of the study area are less prominent with respect to surrounding units, consistent with previous observations suggesting that glaciofluvial formations are often buried under younger sediments in this area (GSC maps; Veillette 1996; Cloutier et al. 2015). The Cadillac Fault, separating the Abitibi (north) and Pontiac (south) geological sub-provinces and roughly corresponding to the position of the water divide between lakes Ojibway (north) and Barlow (south), is used for partitioning the data in Figure 11, allowing consistency with both bedrock characteristics and glacial/post-glacial accumulation processes.

Figure 11. Average unconsolidated sediment thickness (UST, solid line) and glaciofluvial sediment thickness (dashed line) against bedrock surface elevation. The Cadillac Fault, separating the Abitibi (north) and Pontiac (south) geological sub-provinces and roughly corresponding to the position of the water divide between lakes Ojibway (north) and Barlow (south), is used for partitioning the data.

Aquifer delineation: flow conditions and groundwater prospecting

The estimated volumes of confined and unconfined aquifers are reported in Table 5, whereas aquifer thicknesses are shown in Figure 9(d). Overall, assuming that sand and gravel deposits are characterized by a porosity of approximately 30–40% (e.g. see Morris and Johnson 1967; Todd and Mays 2005), it is estimated that the confined and unconfined aquifers of the region store 14–19 km³ and 6–8 km³ of groundwater, respectively. The thicknesses of confined granular aquifers can reach up to 71 m. The greatest aquifers (in terms of thickness and areal extent) seem to be preferentially associated with the northern portion of the study area, where the greatest UST values were also recorded (Figure 9(a)). One regional trend is that the thickest confined aquifers are often found on the flanks of eskers and moraines where glaciofluvial sand/gravel deposits are covered by the regional clayey aquitard (stratigraphic sequences #4, 6, 9 and 13). Similarly, the unconfined aquifers can locally reach 90 m in thickness within the region. The thickest unconfined aquifers generally correspond to the widest eskers and moraines found in the eastern portion of the study area, in the vicinity of Amos and Val-d’Or. From a physical perspective, it can be assumed that the aquifer sectors (both confined and unconfined) characterized by the greatest transmissivities over the widest areas will correspond to the sectors where the groundwater resource is most available in terms of exploitable volumes. The transmissivity (T; [L²T⁻¹]) of an aquifer is defined as the product of its hydraulic conductivity (K; [LT⁻¹]) and its saturated thickness (b; [L]):(8)

Figure 9(d) provides a tool for groundwater prospecting because it allows estimating the b term in Equation (8)(8) for each 100 m × 100 m parcel of the study area, and identifying the widest sectors associated with thick aquifers. The available dataset for documenting the K of unconsolidated sediments at the regional scale (e.g. Cloutier et al. 2016) suggested highly variable values, precluding the construction of a robust transmissivity map based on Figure 9(d). Nevertheless, it provides a tool for targeting areas where field measurements aimed at evaluating K values should be conducted for mapping the availability of groundwater resources.

Discussion

The information obtained from the UST and FGGUT calculations is used jointly with the stratigraphic sequence map to propose an interpretation of the regional hydrogeological framework and to discuss hydrological processes on the basis of unconsolidated units’ architecture. Ultimately, the stratigraphic sequence map is presented as a tool for supporting knowledge transfer to political decision makers and stakeholders in an effort to promote the establishment of strategies oriented toward groundwater protection and sustainable development. Its advantages and limitations are briefly discussed.

The stratigraphic sequence map: A tool for understanding hydrological processes

When working at the regional scale, an understanding of the geological framework is generally central for documenting hydrological processes and the state of the groundwater resource (e.g. see Ross et al. 2005). From recharge to discharge areas, groundwater flows according to hydraulic head gradients within aquifers bounded by aquitards, and a proper representation of the geologic environment is required for interpreting flow processes. In that sense, the stratigraphic sequence map provides a regional-scale representation of the architecture of aquifers and aquitards that can be used to deduce qualitative and quantitative information in relation to groundwater recharge, flow, and discharge along with surface water/groundwater interactions.

The interpretation of the spatial distribution (Figure 9c) and average elevation (Table 5) of the 15 stratigraphic sequences found within the region is used to create a conceptual interpretation of regional flow paths, in relation with the architecture of unconsolidated deposits (Figure 12). The underlying assumption in the use of stratigraphic sequences for delineating recharge/discharge areas is that these processes are most effective within sectors of unconfined aquifers (where the FGGU is absent). Owing to its low hydraulic conductivity, the clay layer is most likely to act as a sub-horizontal barrier limiting infiltration and recharge processes and routing precipitation inputs toward the hydrographic network. Conversely, in sectors where the clay layer is absent, and especially where coarse grained (sand/gravel) sediments are found (stratigraphic sequence #3), water inputs from precipitation are most likely to infiltrate in soils and ultimately become recharge, depending on the losses associated with evapotranspiration. This interpretation is consistent with the pattern of the hydrographic network within the study region, as streams and rivers are rarely found directly on the central portions of eskers and moraines but rather generally find their origin at the margins of eskers/moraines and bedrock/till outcrops, in the vicinity of the unconfined/confined aquifer transition zone (Table 4; stratigraphic sequences #6-9 and 11). This pattern is interpreted herein as reflecting the effect of shallow groundwater recharge–discharge systems that are mainly associated with eskers and moraines. Owing to the coarse-grained sediments composing the eskers and moraines, the water originating from precipitation preferentially infiltrates within soils. Part of this water is lost through evapotranspiration processes, and the remainder contributes to groundwater recharge. Under a permanent regime, it is assumed here that recharge is sufficient for allowing groundwater levels to reach the elevation imposed by the clay aquitard and/or bedrock outcrops on the flanks of eskers and moraines. Most likely in response to a lower transmissivity of confined vs. unconfined aquifers in the region, part of the recharge waters associated with unconfined aquifers (especially eskers and moraines) are not allowed to flow through the clay plain. Under such circumstances, the external margin of the main esker and moraine segments of the region correspond to areas of groundwater discharge feeding extensive peatland complexes and point-specific springs and streams, as conceptualized in Figure 12. Once discharging groundwater is transferred to the hydrographic network, further surface water/groundwater interactions are most likely to be limited over stratigraphic sequences #4–7, 9–10 and 13–14 as the regional clayey aquitard prevents hydrological exchanges. Ultimately, the groundwater recharged in areas laying above the clay plain most likely discharges into streams and rivers in sectors of unconfined aquifers (stratigraphic sequences #1-3, 8, 11-12 and 15) and into the large lakes associated with bedrock depressions. This water subsequently flows toward the James Bay or the St. Lawrence River through rivers.

Figure 12. Conceptual representation of regional distribution of stratigraphic sequences and their position among regional flow paths. The stratigraphic sequences are depicted for illustrating their general distribution and elevation within the regional hydrogeologic environment and the numbers correspond to those found in Tables 4 and 5. The conceptual model suggests increasing unconsolidated sediment thickness (UST) toward the right side of the cross section.

The stratigraphic sequence map: A tool for supporting decision makers

Local and regional political decision makers often must deal with the establishment of land management strategies that lead to a compromise between economic development and environmental protection. Incorporating groundwater protection within this decision scheme is not straightforward because the resource is not visible at the surface and groundwater protection can represent a particular challenge for decision makers (because the resource is not visible at the surface, in opposition to surface waters and ecosystems). The stratigraphic sequence map nevertheless provides a framework for explaining linkages between the geological/geomorphological features of a region and the associated groundwater resources. It could therefore constitute a valuable tool for favoring knowledge transfer to political decision makers and stakeholders dealing with groundwater management and protection.

Protected areas: Ensuring safe access to groundwater for future generations

One central approach for ensuring environmental protection is the identification of protected conservation areas where potentially harmful human activities are prohibited. Such protected areas are often defined in an effort to protect valuable ecosystems and vulnerable species. As a complement, the protection of aquifers could present a unique opportunity to ensure (1) safe access to potable water for future generations, and (2) the maintenance of groundwater functions within ecosystems. In that sense, the stratigraphic sequence map (Figure 9(c)) provides a framework for identifying areas where groundwater protection seems critical. Within the study region, it is proposed that the areas located in the vicinity and on the largest unconfined aquifers found in glaciofluvial complexes (Figure 9(d); Table 4; stratigraphic sequences #3, 8 and 11) should be targeted as potential conservation areas. Three main arguments are proposed to support this statement:

(1)	The eskers and moraines of the region are widely recognized as productive aquifers supplying high-quality groundwater (Nadeau et al. 2015; Cloutier et al. 2016). Protecting the widest of these aquifers, whether they are already exploited or not, would likely stand as a valuable approach for ensuring safe access to potable water for future generations;
(2)	The point-specific springs and diffuse groundwater seepage zones found at the margin of eskers and moraines are associated with the development of groundwater-dependent ecosystems (GDE). One major shortcoming within the Abitibi-Témiscamingue region is that GDEs remain poorly documented and the importance of groundwater functions for surface ecosystems is not yet quantified. Nevertheless, given that such GDEs have been described in other regions where eskers are found (e.g. see Smerdon et al. 2012), it seems fair to suggest that as a precaution, some sectors associated with groundwater discharge areas should be targeted for protection;
(3)	The data, as conceptually represented in Figure 12, reveal that stratigraphic sequences #3, 8 and 11 are generally located in the highest portions of regional flow paths. In that sense, protecting areas associated with these stratigraphic sequences would most likely contribute (at least indirectly) to the protection of confined aquifers found farther along flow paths.

Based on these arguments, it is proposed here that the unconfined aquifers associated with glaciofluvial complexes (stratigraphic sequences #3, 8 and 11) deserve specific protection measures in regional land management strategies. The stratigraphic sequence map (Figure 9(c)) and the conceptual representation of regional flow paths stand as valuable tools for supporting the arguments reported above because they allow representing the hydrogeological environment in an intelligible manner.

Environmental liabilities: Identifying critical issues and defining remediation plans

The socioeconomic context of the Abitibi-Témiscamingue region is largely influenced by the exploitation of natural resources. Forestry and mining stand out as two major regional economic activities. Overall, 794 sand and gravel extraction sites, 57 closed municipal dumps and 127 contaminated sites are identified within the study region (Cloutier et al. 2016). The scale of these environmental liabilities, when compared to the regional population, suggests that applying widespread and synchronous remediation measures at all impacted sites would be unrealistic. Remediation actions should rather be conducted following a stepwise procedure. By way of example, the volumes, ages and nature of contaminants found within contaminated sites along with their distance from private and/or municipal groundwater wells could be taken as criteria for prioritizing remediation actions. In addition, the stratigraphic sequence map (Figure 9(c)) and the regional evaluation of unconsolidated aquifer thicknesses (Figure 9(d)) could provide complementary tools for identifying sectors where remediation should be prioritized. The sectors associated with the thickest shallow unconfined aquifers (Table 5; stratigraphic sequences #3, 8 and 11) should be given special attention because these areas are likely to store significant volumes of groundwater and to be located in the uppermost recharge areas of regional flow paths (Figure 12). In that sense, the stratigraphic sequence map provides a simple tool for decision makers to identify environmental liabilities located in the vicinity of sensitive aquifers.

Advantages and limitations of the stratigraphic sequence map and associated calculations

Advantages of the proposed approaches

The approach used for calculating UST and FGGUT is relatively simple and straightforward, and allows representing USTs realistically on the basis of large datasets. It is especially fit for shield regions characterized by a rugged bedrock surface creating numerous outcrops. Despite its limits, the DNBO functions allow estimating UST between outcrops in areas where borehole data are not available. In that sense, the DNBO functions could stand as a useful complement to interpolation procedures.

The stratigraphic sequence map stands as a hybrid between 2D and 3D representations of the hydrogeological environment. It provides a framework allowing a continuous depiction of the surface geology in the x–y plane coupled with a qualitative representation of the architecture of geological units along the z-axis, without providing information with respect to their thicknesses. The main advantage associated with the stratigraphic sequence map is that it allows visualizing the architecture of buried geological units in a comprehensive and intelligible way at the regional scale. In that sense, it could stand as a highly valuable complement to 3D representations, especially in vast regions (i.e. >10,000 km²) where 2D maps allow a simplified view of the geological environment. Overall, the stratigraphic sequence map provides a communication tool (1) for scientists to share hydrogeological interpretations and (2) for politicians to support decisions related to land management issues. Another significant advantage of the stratigraphic sequence map is that the concept is broadly applicable for representing the architecture of unconsolidated geological units. The visualization concept provided by the map does not depend on the approaches employed for evaluating unconsolidated unit thicknesses and extent. Therefore, different geological models could be represented through the use of a stratigraphic sequence map, notwithstanding the approaches that were used for constructing the model.

Limitations of the approach: Regional data vs. local observations

The different calculation approaches developed in this study for estimating UST and FGGUT (Equations 1 and 5) are specifically based on the exploited extensive regional dataset. In that sense, the approach and equations need to be adapted before being used in other regions. The UST and FGGUT calculation approaches (based on DNBO functions) are best suited for areas of broad glacial lake sediments or till plains that are punctuated by bedrock outcrops. The proposed multiplication factors (MF; Equations 2 and 6) are mainly useful in areas where the bedrock surface presents a significant roughness.

The DNBO functions assume increasing UST with increasing distances from bedrock outcrops. Although the regional dataset supports this assumption, in the absence of nearby borehole data, the UST and FGGUT values remain first-order estimates. DNBO functions are likely to be highly variable in space and will have to be established in other areas. Overall, the DNBO functions do not allow a realistic evaluation of UST for long and narrow features such as valleys and eskers, and corrections need to be applied to obtain adequate estimates for such features (e.g. see Figure 5 – steps 7 and 8). More sophisticated interpolation techniques could improve the evaluation of UST for glaciofluvial formations. However, the lack of borehole data available for characterizing these formations will remain a constraint. Although surface deposits are depicted by the stratigraphic sequences, the specific geological setting may influence the suitability of the stratigraphic sequence map. The approach presented in this study is best suited for an environment where the regional stratigraphy can be simplified to few units. In a context where stratigraphy would be more complex (i.e. numerous layers), the exercise can be expected to be much more complex.

The representation of the geological framework through the use of the stratigraphic sequence map also contains a certain level of uncertainty. First, the use of a 100 m × 100 m mesh size for representing the environment has the advantage of simplifying regional-scale GIS-based calculations, but also leads to significant local-scale inaccuracies. In addition, the creation of the different layers that are required for constructing the map relies on the use of various stratigraphic borehole data, presenting a level of uncertainty that is hard to quantify. Several factors can reduce the reliability of borehole data for constructing a geological model, including (1) compilation errors, (2) inaccuracies in geographic coordinates, and (3) errors in the identification of stratigraphic layers, among others. The aforementioned sources of error can further be magnified through the geoprocessing steps. Such errors are likely to affect the quality of the representations provided by the stratigraphic sequence map. In addition, owing to a lack of data at the regional scale, the maximum thickness of sublittoral sand and peatland was set to 5 m. Although this value is realistic within the regional framework, it represents a source of error within the proposed approach. Developing a quantitative function for estimating the thickness of these deposits would improve the proposed approach. Finally, the spatial distribution of the error associated with the proposed approach was not evaluated here. In that sense, although the approach proposed herein provides a useful tool for interpreting regional-scale hydrological processes, local interpretations should systematically rely on the acquisition of further field data in an effort to corroborate (or invalidate) the information that can be derived from the stratigraphic sequence map.

Conclusions

The general objective of this study was to present an approach allowing the intelligible representation of 3D hydrogeological data as a tool for discussing hydrogeological processes and supporting political decision makers dealing with issues related to groundwater resource assessment, protection and sustainable development. A GIS-based approach was first developed for evaluating the thickness of unconsolidated sediments at the scale of a 19,397 km² region located in Abitibi-Témiscamingue (western Québec, Canada) within the Canadian Shield. These data were subsequently used jointly with surface geology maps in order to construct a stratigraphic sequence map, allowing a continuous depiction of the surface geology in the x–y plane coupled with a qualitative representation of the architecture of geological units along the z-axis. From a scientific standpoint, the stratigraphic sequence map was discussed for:

(1)	Delineating the regional confined and unconfined aquifers;
(2)	Identifying groundwater recharge and discharge areas;
(3)	Constructing a conceptual model of regional groundwater flow paths.

From a social/political standpoint, the stratigraphic sequence map was used as a tool for:

(1)	Identifying areas where groundwater protection should be included in land management constraints in an effort to ensure the safe access to potable water for future generations;
(2)	Providing a framework for prioritizing remediation actions in relation with environmental liabilities.

Overall, the stratigraphic sequence map represents a valuable complement to 3D geological models because it helps to visualize the architecture of buried geological units in a comprehensive and intelligible way at the regional scale. It stands as a valuable tool for facilitating the transfer of hydrogeological knowledge from scientists to political decision makers. Such tools seem critically needed as scientists legitimately tend to focus on the use of models that are growing in complexity and sophistication for improving hydrogeological interpretations.

Acknowledgements

This research was supported by the partners of a regional initiative for sustainable development of groundwater in Abitibi-Témiscamingue, including the Government of Québec (Ministère du Développement durable, de l’Environnement et de la Lutte contre les changements climatiques; Ministère de l’Économie, de l’Innovation et des Exportations; Ministère des Affaires municipales et de l’Occupation du territoire), and the Regional Conference of Elected Officials of Abitibi-Témiscamingue (CRÉAT). We acknowledge the Québec Ministry of Natural Resources (Ministère de l’Énergie et des Ressources naturelles) for their contribution during field work. The authors also acknowledge the financial contribution of the Québec Ministry of the Environment (Ministère du Développement durable, de l’Environnement et de la Lutte contre les changements climatiques) through the Groundwater Knowledge Acquisition Program (PACES program), with significant contributions from regional partners involved in the PACES program, including the Regional County Municipalities (Abitibi, Vallée-de-l’Or, Abitibi-Ouest, Ville de Rouyn-Noranda, Témiscamingue) and the CRÉAT. The authors thank Magalie Roy, Daniel Blanchette, Pierre-Luc Dallaire and Dr. Jean Veillette for their extensive scientific advice for improving the approaches presented herein. The authors kindly thank Dr. Martin Ross and one anonymous reviewer for numerous constructive comments that greatly contributed to improving this manuscript.